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ABSTRACT 


Two  methods  of  estimating  the  attitude  position  of  a  spacecraft  are  examined  in  this 
thesis:  the  extended  Kalman  filter  (EKF)  and  the  unscented  Kalman  filter  (UKF).  In 
particular,  the  UnScented  QUatemion  Estimator  (USQUE)  derived  from  [4]  is 
implemented  into  a  spacecraft  model.  For  generalizations  about  the  each  of  the  filters,  a 
simple  problem  is  initially  solved.  These  solutions  display  typical  characteristics  of  each 
filter  type.  The  UKF  is  very  attractive  in  spacecraft  attitude  estimation,  given  that 
spacecraft  dynamics  are  highly  nonlinear.  For  nonlinear  systems,  the  UKF  is  of 
particular  interest  because  it  uses  a  carefully  selected  set  of  sample  points  that  more 
accurately  map  the  probability  distribution  than  the  linearization  of  the  standard  extended 
Kalman  filter.  This  leads  to  faster  convergence  of  the  attitude  solution  from  largely 
inaccurate  initial  conditions.  The  filter  created  in  this  thesis  is  formulated  based  on 
Markley  and  Crassidis’s  work  on  standard  attitude-vector  measurements  using  a  gyro- 
based  model  for  attitude  propagation.  From  the  standard  attitude  vector  measurements, 
the  global  attitude  parameterization  is  found  and  given  by  a  quaternion,  while  a 
generalized  three-dimensional  attitude  representation  is  used  to  define  the  local  attitude 
error.  The  multiplicative  quaternion-error  is  then  found  from  the  local  error.  The 
simulation  results  indicate  that  the  unscented  filter  is  more  robust  than  the  extended 
Kalman  filter. 
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I.  INTRODUCTION 


A.  BACKGROUND 

In  August  of  1960,  the  United  States  Air  Force  and  the  Central  Intelligence 
Agency  successfully  launched  the  world’s  first  reconnaissance  satellite,  Corona.  The 
imaging  resolution  was  8  meters  and  taken  on  film.  The  program  lasted  for  12  years,  and 
ushered  in  the  era  of  space-based  reconnaissance  and  intelligence  gathering  that  would  be 
iconic  of  the  Cold  War.  Since  the  beginning  of  spacecraft  building,  organizations  have 
prized  themselves  on  pushing  the  envelopes  of  technology.  Programs  such  as  NASA’s 
Explorer,  TIROS,  and  Pioneer  later  proved  that  the  U.S.  was  investing  heavily  on  space 
technologies. 

I 


Figure  1.  Dr.  William  Pickering,  Dr.  James  Van  Allen,  and  Dr.  Wemher  Von  Braun 
hold  a  model  of  the  Explorer  1  vehicle  above  their  heads.  Credit:  NASA 

Historically,  the  building  of  spacecraft  has  been  a  lengthy  process,  often  taking 
many  years  or  even  decades.  Recently,  however,  a  new  methodology  for  building 
spacecraft  has  transpired.  Organizations  such  as  The  Office  of  Responsive  Space  (ORS) 
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have  been  created  to  change  this  expensive  and  lengthy  process  into  one  that  focuses  on 
providing  a  "good  enough"  service  in  a  timely  manner  [1],  This  push  for  faster  programs 
has  also  led  areas  of  the  industry  to  build  smaller  systems  in  attempts  to  utilize  the 
leftover-over  payload  mass  of  launch  vehicles.  This  is  more  commonly  referred  to  as  a 
"secondary  payload."  The  industry  push  combined  with  the  emerging  university 
nanosatellite  community  has  created  an  influx  of  new  commercialism  for  space-based 
hardware. 

One  of  the  limiting  technologies  in  the  small  spacecraft  arena  is  attitude 
determination  and  control  systems  (ADCS).  While  currently  there  is  an  increased  interest 
in  this  area,  a  limited  number  of  complete  solutions  in  a  3U  or  lU-class  nanosatellite 
have  been  demonstrated  on-orbit.  Many  proposed  solutions  are  also  not  affordable  to  this 
community.  While  companies  like  Boeing,  Honeywell,  and  Sinclair  are  working  on 
hardware  solutions,  the  problem  of  attitude  determination  and  control  can  be  attacked 
from  both  sides.  That  is  to  say,  as  the  hardware  is  being  developed,  both  academic  and 
commercial  institutions  can  focus  their  resources  on  the  optimal  estimation  and  control 
theory  problems.  The  lack  of  an  affordable  hardware  should  not  inhibit  willing  parties  to 
develop  solutions  and  methods  for  the  small  spacecraft  ADCS  problem. 

B.  SHORT  OVERVIEW  OF  ATTITUDE  ESTIMATION 

One  of  the  most  common  estimation  techniques  that  has  been  widely  used  for 
various  dynamics  systems  is  the  Kalman  filter.  While  the  filter  was  initially  designed  for 
linear  systems,  variations  of  this  filter  have  been  developed  in  particular  for  nonlinear 
systems.  The  extended  Kalman  filter  can  be  used  on  nonlinear  systems  and  is  based  on 
linearizing  the  system  dynamics.  While  this  is  potentially  attractive  for  the  nonlinear 
spacecraft  attitude  control  problem  there  are  several  associated  with  the  nonlinearization 
[2]. 

While  the  EKF  has  proven  to  be  a  popular  tool  for  nonlinear  estimation,  it 
continues  to  endure  some  fundamental  issues  inherent  in  the  linearization  process,  which 
can  be  the  potential  cause  of  divergence.  A  later  development  for  nonlinear  estimation 
was  developed  by  Julier  and  Uhlmann  and  is  called  the  “unscented”  Kalman  filter  [3]. 
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The  UKF  is  “founded  on  the  intuition  that  it  is  easier  to  approximate  a  probability 
distribution  than  it  is  to  approximate  an  arbitrary  nonlinear  function  or  transformation 
[3]."  The  UKF  successfully  avoids  the  EKF  linearization  step  by  introducing  a  set  of 
sample  points  that  capture  the  higher  order  statics  of  the  system.  Finally,  the  UKF 
method  has  been  developed  to  estimate  the  quaternions  associated  with  the  attitude  of  a 
spacecraft  [4],  The  numerical  simulations  presented  in  these  studies  have  illustrated  the 
superior  performance  of  the  UKF  in  this  context. 

The  primary  goal  of  this  thesis  is  to  develop  and  verify  estimation  algorithms  and 
simulation  code  for  a  spacecraft  attitude  determination  system  (ADS).  In  particular,  the 
two  estimation  methods  that  are  compared  for  determining  the  attitude  are  the  extended 
Kalman  fdter  (EKF)  and  the  unscented  Kalman  filter  (UKF).  Each  fdter  is  evaluated 
based  on  error  computation  time.  The  inherent  linearity  and  nonlinearity  of  each  type  of 
fdter  is  examined  by  choosing  related  problems  that  highlight  issues  in  trying  to  use  a 
linear  fdter  (EKF)  to  solve  a  nonlinear  problem.  To  do  this,  two  separate  simulations 
codes  were  designed.  These  simulation  codes  include  an  accurate  spacecraft  model  where 
torque  disturbances,  Earth  physics,  and  orbital  mechanics  are  accounted  for,  as  well  as 
sensor  models  of  an  inertial  measurement  unit  and  magnetometer. 

A  simplified  problem  was  used  to  verify  the  behavior  of  both  estimation  methods 
on  linear  and  nonlinear  dynamics.  For  this,  the  simple  pendulum  was  used  as  a  way  to 
show  how  each  fdter  can  be  used  to  estimate  the  states  of  a  given  dynamic  problem. 
After  this  problem  was  worked,  these  filters  were  used  as  analogs  against  a  simulated 
spacecraft  model.  Primarily  focusing  on  the  UKF,  this  thesis  discusses  the  differences 
between  the  two  filters  and  focuses  on  the  benefits  of  using  nonlinear  estimation.  It  is 
widely  known  that  there  are  many  benefits  to  nonlinear  estimation.  The  UKF  is  very 
attractive  in  spacecraft  attitude  estimation,  given  that  spacecraft  dynamics  are  highly 
nonlinear.  This  thesis  highlights  these  benefits  while  solving  both  the  EKF  and  UKF 
spacecraft  attitude  estimation  problem.  While  previous  theses  discussed  the  nuances  of 
characterizing  these  types  of  sensors  for  inclusion  in  the  simulation  [4],  this  paper  will 
focus  on  the  estimation  methods  as  they  apply  to  attitude  determination. 
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C.  RECENT  CUBESAT  ADCS  SYSTEMS 

Several  spacecraft  that  have  implemented  ADCS  systems  into  the  small  CubeSat 
standard.  This  section  discusses  three  of  these,  which  include  the  following: 

•  Canadian  Advanced  Nanospace  experiment  (CanX)  -The  University  of 
Toronto  Institute  for  Aerospace  Studies  Space  Flight  Laboratory  (UTIAS 
SLF) 

•  AISSat-1  -  The  University  of  Toronto  Institute  for  Aerospace  Studies 
Space  Flight  Laboratory  (UTIAS  SLF) 

•  Radio  Aurora  experiment  (RAX)  -  The  University  of  Michigan 

1.  Canadian  Advanced  Nanospace  experiment  (CanX)  ADS 

The  Canadian  Advanced  Nanospace  experiment  (CanX)  program  is  run  by  the 
University  of  Toronto  Institute  for  Aerospace  Studies  (UTIAS)  Space  Flight  Laboratory. 
CanX-1  launched  on  June  2003  from  Plesetsk,  Russia,  and  was  a  1-U  CubeSat  that 
consisted  of  several  ADS  hardware  components.  The  primary  mission  of  CanX-1  was  to 
demonstrate  the  several  experimental  ADS  components.  The  CanX-1  ADS  package 
consisted  of  a  CMOS  Imager  for  ground-controlled  horizon  sensing  and  star  tracking, 
active  three-axis  magnetic  stabilization  and  a  Global  Positioning  System  (GPS)  receiver 
that  was  modified  to  work  in  low  Earth  orbit.  Figure  2  shows  a  picture  of  the  CanX 
CMOS  imager  used  for  star  tracking  [6], 
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Figure  2. 


CanX-1  Agilent  Technologies  CMOS  Imager 


CanX-2,  which  launched  in  April  2008,  uses  many  of  the  same  types  of  ADS 
systems.  The  CanX-2  ADS  uses  a  suite  of  sun  sensors  and  a  three-axis  magnetometer. 
Both  CanX-1  and  CanX-2  use  a  standard  extended  Kalman  filter  to  estimate  the  attitude 
of  the  spacecraft. 


Figure  3.  Computer  Rendering  of  CanX-2 
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2. 


AISSat-1  ADS 


AISSat-1  is  a  6-kg  Norwegian  nanosatellite,  being  constructed  on  behalf  of 
government  of  Norway  by  UTIAS/SFL,  whose  primary  mission  is  to  investigate  the 
feasibility  and  performance  of  a  spacecraft-based  Automatic  Identification  System  (AIS) 
sensor  in  low-Earth  orbit  as  a  means  of  tracking  maritime  assets.  AISSat-1  is  intended  as 
both  a  research  and  development  platform,  and  a  demonstration  mission  for  a  larger 
operational  capability. 


Figure  4.  Computer  Rendering  of  AISSat-1,  from  [7] 


A  full  3-axis  attitude  determination  and  control  system  provides  attitude 
stabilization  and  fine  pointing  for  AISSat-1.  The  satellite  is  able  to  point  in  either  and 
inertial  orientations,  or  an  orbit-frame-fixed  orientation,  including  on  nadir.  Attitude 
sensors  consist  of  six  sun  sensors,  a  magnetometer  and  rate  gyros.  Three  orthogonally 
mounted  reaction  wheels  and  three  magnetorquer  coils  controls  the  actuation  of  the 
satellite.  The  magnetorquer  is  used  for  de-tumbling  and  momentum  dumping  while  the 
reaction  wheels  provide  fine  attitude  pointing  capability.  The  attitude  control  system  is 
able  to  maintain  several  degree  level  pointing  accuracy  and  stability  over  the  course  of 
the  entire  orbit,  including  eclipse.  For  attitude  estimation,  this  spacecraft  also 
implemented  an  extended  Kalman  filter  [7]. 
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3.  Radio  Aurora  Explorer  (RAX)  ADS 

The  Radio  Aurora  Explorer  (RAX)  spacecraft,  currently  being  developed  by  The 
University  of  Michigan  is  a  3U  CubeSat,  which  will  also  implement  and  attitude 
determination  system.  The  primary  scientific  objective  of  the  Radio  Aurora  Explorer 
(RAX)  mission  is  to  understand  the  microphysics  of  plasma  instabilities  that  lead  to  field- 
aligned  irregularities  (FAI)  of  electron  density  in  the  polar  lower  (80-400  km) 
ionosphere.  For  attitude  control,  an  inertial  measurement  unit  in  conjunction  with  sun 
sensors  and  magnetometers  will  observe  the  time  it  takes  the  passive  magnetic  attitude 
control  system  to  de -tumble  the  spacecraft  after  deployment.  This  system  will  implement 
a  continuous-discrete  extended  Kalman  filter.  They  will  implement  a  13  state  filter, 
which  will  consist  of  3  position,  3  velocity,  4  quaternions,  and  3  angular  rates.  The  team 
will  implement  the  QUatemion  ESTimator  (QUEST)  method  developed  by  Shuster  and 
Oh  [8].  Literature  describes  the  QUEST  method  as  computationally  expensive;  however, 
the  information  will  be  gathered  on  orbit  and  processed  on  the  ground  to  eliminate 
computational  constraints  on  the  filtering  process.  Some  of  the  ADS  hardware  will 
include  six  3 -axis  magnetometers,  nine  sun  sensors  and  an  inertial  measurement  unit, 
which  will  consist  of  a  3 -axis  gyroscope. 
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II.  GENERAL  EXTENDED  AND  UNSCENTED  KALMAN 
FILTERING  METHODS  FOR  THE  ESTIMATION  OF 
DYNAMIC  SYSTEMS 


A.  BACKGROUND 

Accurate  attitude  knowledge  is  essential  for  many  spacecraft  missions.  Kalman 
fdtering  has  been  widely  known  since  the  1960s  as  a  method  for  filtering  out  noise  in  a 
given  measurement.  Theoretically,  the  Kalman  filter  is  a  sequential  optimal  estimator  for 
what  is  called  the  linear-quadratic  problem,  which  is  the  problem  of  estimating  the 
instantaneous  "state"  of  a  linear  dynamic  system  including  its  uncertainty— by  using 
measurements  linearly  related  to  the  state  corrupted  by  white  noise  [10].  For  attitude 
determination,  several  types  of  Kalman  filters  have  been  developed  over  the  years.  This 
section  describes  two  basic  types  of  Kalman  filters,  the  extended  Kalman  filter  (EKF)  and 
the  unscented  Kalman  filter  (UKF).  There  have  been  many  technical  papers  written  on 
Kalman  filtering  for  state  estimation  [2]  [3].  This  chapter  will  start  the  discussion  with 
the  continuous-time  Kalman  filter  as  a  base  line.  Several  textbooks  use  a  variety  of 
nomenclature  to  describe  this  estimation  process.  As  a  standard,  the  following  tables 
define  the  notation  used  in  this  thesis.  These  tables  are  also  consistent  with  [9]  and  [10], 
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Table  1. 


Standard  Symbols  of  Kalman  Filtering,  from  [10] 


Symbols 

Symbol  Definition 

F 

Dynamic  coefficient  matrix  (state  matrix)  of  a  continuous  linear 

differential  equation  defining  a  dynamic  system 

G 

Coupling  matrix  between  random  process  noise  and  the  state  of  a  dynamic 

system 

H 

Measurement  sensitivity  matrix  defining  the  linear  relationship  between 

the  state  of  the  dynamic  system  and  measurements  that  can  be  made,  (also 

known  as  a  coefficient  matrix  [9]) 

K 

Kalman  gain  matrix 

P 

Covariance  matrix  of  state  estimation  uncertainty 

Q 

Covariance  matrix  of  process  noise  in  the  system  state  dynamics  also 

called  the  process  noise  covariance 

R 

Covariance  matrix  of  observational  (measurement)  uncertainty  also  called 

the  measurement  noise  covariance 

X 

State  vector 

y 

Vector  (or  scalar)  of  measured  values. 

<D 

State  transition  matrix  of  a  discrete  dynamic  system 
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Table  2.  Special  State  Space  Notation,  from  [10] 


Symbols 

Symbol  Definition 

xkd) 

The  z-th  component  of  the  vector  x,  or  the  z-th  element  of  the  sequence.  The 
sub-index  k  refers  to  the  sequence  of  propagation  as  it  occurs  in  the  filtering 
process,  i.e.  k+1  can  be  referred  to  as  the  “update”  term  that  is  determined 
from  the  same  term  calculated  previously.  xk+1  ( i )  =  xk  (i i )  +  noise 

X 

An  estimate  of  the  value  of  x. 

K 

A  priori  estimate  of  the  xk ,  conditioned  on  all  prior  measurements  except 
the  one  at  time  tk 

X  + 
xk 

A  posteriori  estimate  of  the  x ,  conditioned  on  all  available  measurements  at 
time  tk 

y 

A  measurement  of  some  quantity  we  can  estimate  to  the  state  vector  from. 

X 

Derivative  of  x  with  respect  to  time 

The  Kalman  filter  uses  a  parametric  characterization  of  the  probability 
distribution  of  its  estimation  errors  in  determining  the  optimal  filtering  gains,  and  it  is  the 
probability  distribution  that  can  be  used  for  assessing  its  performance  as  a  function  of  the 
“design  parameters”  of  an  estimation  system  [10].  Some  of  these  can  include: 

•  the  types  of  sensors  used, 

•  the  locations  and  orientations  of  the  various  sensor  types  with  respect  to 
the  system  to  be  estimated, 

•  the  allowable  noise  characteristics  of  the  sensors, 

•  the  data  sampling  rates  for  the  various  sensor  types,  and  most  importantly, 

•  the  level  of  model  simplification  to  reduce  implementation  requirements. 
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B.  CONTINUOUS-TIME  EXTENDED  KALMAN  FILTER 

For  nonlinear  systems,  such  as  spacecraft  dynamics,  the  extended  Kalman  filter 
(EKF)  has  been  previously  proposed  in  literature  and  used  on-board  many  spacecraft  [2], 
In  the  EKF,  the  state  transition  and  observation  models  do  not  need  to  be  linear  functions 
representing  the  state,  granted  they  are  differentiable.  Given  that  a  vast  majority  of 
nonlinear  problems  can  be  described  with  differentiable  nonlinear  functions,  the 
Continuous-Time  EKF  can  often  be  used.  The  Continuous-Time  EKF  is  very  similar  to 
the  Continuous-Time  Linear  Kalman  filter  [9].  The  derivation  of  the  Continuous-Time 
EKF  starts  with  the  continuous  non-linear  system  model  below: 

x(t)  =  f(x(t),u(t),t)  +  G(t)w(t)  2.1 

y(0  =  h(x(t),t)  +  v(0  2 .2 

where  it  is  important  to  note  that  f(x(t),  u(t),t)  represents  nonlinear  continuous  function  or 
the  state  transition  model,  while  G(t)  and  w(t)  represent  the  coupling  matrix  and 
continuous-time  covariance  respectively.  For  Equation  2.2,  y (t)  represents  the  measured 
nonlinear  observed  model  using  a  continuous  function  h(x(t),t)  plus  the  continuous-time 
covariance,  v(t). 

The  inherent  linearization  process  can  cause  the  filter  to  diverge,  as  the  Gaussian 
input  does  not  necessarily  produce  a  Gaussian  output  [9],  To  continue,  we  must  assume 
that,  for  our  purposes,  a  linear  representation  of  our  non-linear  system  will  suffice.  For 
example,  this  method  can  certainly  be  used  for  functions  where  small  angle 
approximation  is  valid.  Examples  of  the  limitations  are  discussed  in  detail  in  Section  F. 
For  the  EKF,  we  must  also  assume  that  the  true  state  of  the  system  is  sufficiently  close  to 
the  estimated  state.  Therefore,  the  error  dynamics  can  be  reasonably  approximated  by  a 
linearized  first  order  Taylor  series  expansion.  The  first  order  expansion  of  f(x(t),u(t),t) 
about  a  nominal  state  x(t)  becomes: 
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f  (x(0,  u(0,  t)  =  f  (x(0,  u  (0, 0 + 77 

ox 


[x(0-x(0] 


x(0 


2.3 


where  x(t)  is  close  to  x(t) .  Similarly,  the  output  in  Equation  2.3  becomes  [9]: 


h(x(0,0  =  h(x(0,0  + 


8h 

dx 


[x(0-x(0] 

x(0 


2.4 


Here  the  EKF  solves  this  problem  by  calculating  the  Jacobians  of  /  and  h  around  the 
estimated  state,  which  in  turn  yields  a  trajectory  model  function  centered  around  this 
state.  Figure  5  shows  this  graphically  [11]. 


Figure  5.  Illustration  of  Extended  Kalman  filter  linearization  of  nonlinear  function 

and  the  related  Gaussian  distribution. 

To  find  the  estimate  of  the  state,  the  extended  Kalman  filter  continues  with 
assumption  made  earlier,  that  x(t)  =  x(t) .  Thus,  the  expectation  of  both  Equations  2.3 
and  2.4  gives  the  following  equation,  where  E  represents  the  conditional  mean  or 
expectation  [9]. 
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E  {f  (x(0,  u(0, 0}  =  f <i(0»  u(0, 0 
^{h(x(0,0}=h(x(0,0 


2.5 


2.6 


Therefore,  the  extended  Kalman  filter  for  the  state  and  output  estimate  is  given  by  the 
following  two  equations  [9], 


x(t)  =  f  (x(t),  u(0, 0  +  K(t)[y{t)  -  h(x(t),  0] 
y(0  =  h(x(t),0 


Because  the  equation  of  the  measurement  of  the  state  vector  has  the  same  structure  as  the 
linear  Kalman  filter,  we  can  use  the  covariance  expression  shown  in  Table  3.  The 
following  table  summarizes  the  equations  for  the  continuous-time  extended  Kalman 
filter. 
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Table  3. 


Continuous-Time  Extended  Kalman  Filter,  from  [9] 


Model 

x(0  =  f  (x(0,  u(t),  t )  +  G(t)w(0, 
y(0=h(x(0,0+v(0 

Initialize 

x(c)  =  x0 

P0=E{x(to)xr(t0)} 

Gain 

K{t)  =  P{t)HT(x(t),t)R-\t) 

Covariance 

P(t)  =  F{x(t),t)P(t)  +  P(t)FT(x(t),t ) 
-P{t)HT{x{t),t)R-\t)H{x{t),t)P{t)  +  G(t)Q(t)GT(t) 

^  A  df  x  Sh 

(x,9=s  ■ 

x(0  x(0 

Estimate 

x(t)  =  f  (x(0,  u(0, 0  +  K(t)  [y  (0  -  h(x(0,  o] 

C.  DISCRETE-TIME  LINEAR  AND  EXTENDED  KALMAN  FILTERS 

While  understanding  the  basics  of  the  continuous-time  extended  Kalman  filter  is 
valuable  in  the  sense  that  it  can  often  be  used  to  solve  entire  solutions  analytically, 
implementation  of  this  is  not  practical.  In  most  cases,  the  control  system  is  responding  to 
different  given  inputs.  The  use  of  real-time  processing  is  inevitable  in  the  practical 
implementation  of  estimating  dynamic  systems.  Thus,  the  continuous-time  Kalman  filter 
must  be  discretized  so  that  it  may  be  applied  to  iterative  methods.  This  section  describes 
how  the  Kalman  filter  is  derived. 

Derivation  of  the  discrete-time  filter  and  the  extended  Kalman  filter  are  very 
similar.  To  derive  the  discrete-time  Kalman  filter,  an  assumption  must  be  made  that  both 
the  model  and  measurement  are  available  in  discrete  form.  Here,  we  can  start  with  the 
non-linear  "truth"  model  shown  below  [9]: 


2.10 


y  k=HkM+vk 


where  ^  is  the  state  transition  matrix,  T  is  the  control-input  matrix  that  is  applied  to  the 
control  vector  uk,  and  Y  is  the  noise  matrix.  The  definition  of  ® ^ ,  and  Y  are  shown 
below. 


O  =  eFAl 


2.11 


2.12 


eldt 


2.13 


where  B  and  G  are  the  coefficient  matrices  taken  from  the  continuous  system.  Again,  in 
Equations  2.9  and  2.10  Wk(t)  and  Vk(t)  are  assumed  to  be  zero-mean  Gaussian  white-noise 
processes  and  their  covariance's  are  given  by  the  expectation  equations  [9]: 


£{w,wj} 


k  =  j 


k*j 
k  =  j 


2.14 


2.15 


The  Qk  matrix  accounts  for  the  state  process  noise  while  the  Rk  matrix  accounts 
for  the  expected  measurement  noise.  These  equations  imply  that  the  errors  are  not 
correlated  forward  or  backward  in  time.  We  can  also  assume  that  vk  and  wk  are 
uncorrelated  so: 

£{VX}  =  °  2.16 

Updating  the  current  estimate  of  the  state  xk  to  obtain  xk+]  based  upon  all  k+1 

measurement  subsets  assumes  that  the  gain  (K)  can  vary  in  time.  This  propagation  can 
be  done  using  [9]: 
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2.17 


*M=®X+TlPk 


Furthermore,  the  updated  state  is  given  by: 


*\=K+Kk\yk~H&\ 


2.18 


where  yk  is  the  measurement  vector.  The  gain  Kj(  changes  with  time  properly  weighting 

the  relative  confidence  of  the  accuracy  of  the  propagated  state  verses  the  measured  state. 
To  find  Kk,  first  the  state  error  and  error  covariance  matrixes  must  be  defined  [9] : 


where 


pk  =E{rkrkT] 


Xk=Xk-*k 


2.19 


2.20 


Substituting  Equations  2.9  and  2.17  into  Equation  2.20  and  substituting  the  resulting 
equation  into  Equation  2.19  leads  to: 


P^=<bkK<  +  EkQkYTk 


2.21 


Because  wk  and  are  uncorrelated  the  terms  £'|w/[x^'/  j  =  =  0 .  To  find  the 

updated  error  covariance  matrix,  we  can  use  Equations  2.10  and  2.18.  Then  substitution 
of  the  resulting  equation  into  Equation  2.20  leads  to: 


2.22 


To  find  the  gain  K,  the  trace  of  error  covariance  matrix  Pk  is  minimized.  Solving  gives: 


Kt  =^Hl{rk)[Hk(ri)P^HTl{i-l)+Rl 


1-1 


2.23 
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As  mentioned  previously,  the  extended  Kalman  filter  and  discrete-time  Kalman 
filter  are  nearly  identical.  The  only  difference  between  these  two  are  the  initial  model 
equations  and  the  propagation  equations.  The  extended  Kalman  filter  assumes  that  the 
model  is  a  continuous  function  and  thus  be  differentiable.  This  is  clearly  evident  in  Table 
4. 

Table  4.  Discrete-Time  Linear  Kalman  Filter,  from  [9] 


Model 

x*+i=®*x*+r*u*+Tw*>  wt~AT(o  ,Qk) 

+  W  \k~N(0,Rk) 

Initialize 

2 

<*° 

A  2 

cS  jiL 

11 

a? 

Gain 

Kk  =  Pk~Hl (i- )[Hk  (rk)Pk~HTk  (x- )  +  Rt ]"‘ 

Update 

**=  a +Kt[y* -#*(**)] 

p;=[i-KkHk]p~ 

Propagation 

K+i  =  ®kK+TiPk 

D.  UNSCENTED  KALMAN  FILTER 

The  inherent  issue  with  propagating  Gaussian  random  variables  through  a 
nonlinear  function  can  also  be  approached  using  a  technique  described  as  the  unscented 
transform.  While  the  extended  Kalman  filter  has  many  applications,  and  is  the  most 
popular  method  for  nonlinear  estimation  to  date,  the  unscented  Kalman  filter  (UKF)  was 
proposed  by  Julier,  Uhlmann,  and  Durrant- Whyte  [12]  to  overcome  the  instabilities 
associated  with  the  EKF.  While  the  EKF  typically  works  well  in  the  regions  where  the 
first-order  Taylor  series  linearization  adequately  approximates  the  nonlinear  probably 
distribution,  a  primary  area  of  concern  is  during  the  initialization  stage,  where  the 
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estimated  initial  state  can  be  far  from  the  true  state  [9].  The  UKF  typically  involves  more 
complex  computations  than  the  EKF,  but  has  the  following  advantages: 

1 .  the  expected  error  is  lower  than  the  EKF 

2.  it  can  be  applied  to  non-differentiable  function 

3.  it  avoids  the  derivation  of  Jacobian  matrices 

4.  it  is  valid  to  higher-order  expansions  than  the  standard  EKF  [4] 

The  UKF  can  be  thought  of  as  an  extension  of  the  traditional  Kalman  filter  for  the 
estimation  of  nonlinear  systems  that  implements  the  unscented  transformation.  The 
unscented  transformation  uses  a  set  of  sample,  or  sigma,  points  that  are  determined  from 
the  a  priori  mean  and  covariance  of  the  state.  The  sigma  points  undergo  the  nonlinear 
transformation.  Then  the  a  posteriori  mean  and  covariance  of  the  state  are  determined 
from  the  transformed  sigma  points.  This  approach  gives  the  UKF  better  convergence 
characteristics  and  greater  accuracy  than  the  EKF  for  nonlinear  systems  [13].  The  ability 
of  the  UKF  to  accurately  estimate  nonlinearities  make  it  attractive  for  implementation  on 
spacecraft  as  the  state  and  observations  are  inherently  nonlinear.  This  section  describes 
the  basic  derivation  of  the  unscented  Kalman  filter,  while  the  subsequent  sections 
describe  the  implementation  of  the  UKF  for  attitude  determination. 

The  derivation  of  the  unscented  Kalman  filter  starts  by  selecting  a  nonlinear 
system  defined  by  [9]: 

xt+1=f(x*^)  +  <^w*  2.24 

&=  *>(**,£)  +  ▼*  2.25 

where  Xk  is  the  n  x  1  state  vector  and  yk  is  the  m  x  1  measurement  vector.  It  is 
interesting  to  note  that  a  continuous-time  model  can  also  be  expressed  in  the  form  of 
Equation  2.24.  Similar  to  the  previous  derivations,  Vk  represents  the  measurement-error 
noise  while  Wk  describes  the  white  Gaussian  process  noise  with  covariances  given  by 

E'|w(t)wr(r)|  =  Q(t)5(t  -  t)  2.26 
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E  jw(7)vr(r)}  =  R(t)S(t-r ) 
E  jv(7)wr(r)j  =  0 


2.27 


2.28 


The  covariance  matrices  of  each  of  these  are  given  by  Qk  and  Rk  respectively  [4],  The 
Kalman  filter  update  equations  are  rewritten  from  Table  4  as  [3]: 

K  =  K  +  Kk'>k  229 

K=Pk~-KkPkvvKTk  2.30 


where  vk  is  the  innovations  process,  given  by 

ok  =  yt-y~t=yk-HK>ukk) 


The  covariance  of  the  innovations  process,  vk  is  given  by  Pkvv  [4], 


pvv  —  pyy  D 

rk+ 1  —  rk+ 1  "r  ^£+1 


The  Kalman  gain  is  computed  by  the  following  equation  [4], 


2.32 


Kk  =  Pky(pr)A  2.33 

where  Pf  is  the  cross-correlation  matrix  between  x“,  and  y~ .  The  cross-correlation  is 

defined  later  in  the  discussion  below.  To  define  the  propagation  equations,  the  following 
sigma  points  must  be  computed  [4],  The  filter  starts  by  augmenting  the  state  vector  to  L 
dimensions  in  the  original  state-vector,  model  noise,  and  measurement  noise  where  L  is 
the  size  of  the  vector  x“,  or  the  augmented  state  defined  by  Equation  2.37  [9].  The 

covariance  matrix  is  similarly  augmented  and  this  forms  the  augmented  state  estimate 
vector  shown  below. 
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2.34 


a  <—  2L  columns  from 


z:(  o) 


=  xt 


%°t(i)=<Tk(i)+n 


where  x“  is  an  augmented  state  defined  by  [4] 


2.35 

2.36 


~xk~ 

"  *k  " 

II 

wk 

\a  — 

0^x1 

3_ 

_^mxl  _ 

Augmenting  the  covariance  requires  the  computation  of  2(q+l)  additional  sigma  points. 
It  is  important  to  mention  here  that  q  is  the  dimension  of  wk ,  l  is  the  dimension  of  vk , 
and  m  is  the  output  dimension.  While  L  is  the  size  of  the  vector  x“ ,  the  parameter  y  is 
given  by  the  following  [4], 


y  =  a Jl  +  X  2.38 

and  the  composite  scaling  parameter ,  A,,  is  given  by 

A  =  a\L  +  /c)-L  2.39 

The  constant  a,  represents  the  spread  of  sigma  points  and  is  usually  set  to  a  small  positive 
value  (e.g.,  1x10  4  <  a  <  1).  There  are  2L  values  for  <Jk,  each  representing  the 

positive  and  negative  values  of  the  square  root.  The  Cholesky  method  is  often  used  to 
find  the  square  root  of  a  matrix.  Similar  to  the  EKF,  the  UKF  now  propagates  these 
sigma-points  from  a  Gaussian  distribution  through  a  nonlinear  function,  and  recreates  a 
Gaussian  distribution  by  calculating  the  mean  and  covariance  of  these  results  [11], 
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y=f(x) 


Figure  6.  Illustration  of  the  unscented  Kalman  filter  sigma-points  propagation 


These  sigma  points  are  evaluated  by: 


zk+ 1(0  =  *(xZ(0,x?  (0 


2.40 


where  ^  (/')  is  a  vector  of  the  first  n  elements  of  xl  (0  >  and  Xi'(0  is  a  vector  of  the  next  q 
elements  of  xl  (0 »  with 


Xk  = 


xi  (0 

XkC 0 

Xk(0 


2.41 


The  predicted  mean  for  the  state  estimate  is  calculated  using  a  weighted  sum  of  points 
XL i(0,  given  by: 


*;.,=! y—xiM 


i=  0 


2.42 


where  the  weight  terms  W"ie“"  is  given  by: 
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and 


L  +  A 


2.43 


wmean  -  ■ 

rr  o 


j^rmean  _ yyz°y  _ 


/l 


2(L  +  A) ' 


i  =  l,2,...,2L 


2.44 


Similarly,  the  predicted  covariance  term  is  given  by: 


pm  =  Z  [xZ  (0  -  J_xZfi)  ~  *Z  ]' 

i= 0 


2.45 


where  the  weight  terms  are  given  by  2.44,  and  the  following  equation. 

wr  =^^  +  (\-a2  +  P) 

0  L  +  A  y 


The  mean  observation  is  given  by 


9Z  =Zr*”n.,(0 


;•= 0 


2.46 


2.47 


where 


Y*+i(0  =  h{zZ+i(f),vk+i,ZZ+i(Q>k  +  l) 


2.48 


The  output  covariance  matrix  is  given  by: 


^  =Z»T”[rM(0»<yw][rM(0^yM]r  2,49 

1=0 

The  innovations  covariance  is  given  by  Equation  2.32.  The  cross  correlation  matrix  is 
finally  described  as 


PZ  =  Z  wr  [xZ  (0  -  ]  [n„  (0  -  y:.,  ]' 

1=0 


2.50 
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Finally,  the  Kalman  gain  and  states  are  updated  using  the  following  equations. 

Kk=p-xprr 


pk+=p--KkprKTk 

y*=y*-y*=y*-h(x;,u  kk) 


2.51 

2.52 

2.53 

2.54 


A  summary  of  these  equations  are  listed  in  Table  5  and  will  be  referred  to  in  subsequent 
sections  that  describe  the  implementation  of  these  filters. 
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Table  5. 


Unscented  Kalman  Filter,  from  [9] 


Model 

At+i 

y*  =  *»(**,«*>  vt,*) 

Initialize 

p(K)  =  p0 

Gain 

Kk =pkxy(prr 

Update 

X*  =  x“  +  Kkvk 

p;=pk~-KkprzTk 

^=yk-y~k=yk-HK’ukk) 

Propagation 

i= 0 

p;„ = S  wr  [x’an  -  ]  [xi«  <o  -  J 

i= 0 

z=0 

pz  =  Z, wr  (0  -  fc,  ]  [>„,  (o  - -  y;.,  ]r 

z=0 

_  pTT 

^  £+1  ~~  *  k+\ 

C,  =  Z  wr  [xL  (0  -  ]  [n.,  <0  -  yM  7 

*=0 

E.  IMPLEMENTATION  OF  EKF  AND  UKF  METHODS  USING  THE 
SIMPLE  PENDULUM  PROBLEM 

Prior  to  implementing  the  EKF  and  UKF  on  the  spacecraft  model,  an  easier 
problem  was  solved.  For  this,  the  simple  pendulum  was  used.  Figure  7  shows  a  diagram 
of  the  simple  pendulum  problem. 
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p 


Figure  7.  Simple  Pendulum  Problem 
The  dynamic  equation  is  commonly  known  and  listed  below. 


q  _  ~mgl  sin(<9) 
/,  +  ml2 


Translated  into  the  state  space  model,  this  becomes: 


2.55 


2.57 
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We  can  now  implement  the  state-space  model  into  the  simulation  block  diagram  as  our 
dynamics  state  function.  The  three  states  that  were  estimated  were  0,  0 ,  and  /  . 


Figure  8.  Simulink  Block  Diagram  of  Simple  Pendulum  Model 


By  solving  Equation  2.57,  we  can  then  use  its  solution  to  determine  our 
measurement  equation. 


ymeas=B OSin(0  +  tf) 


2.58 
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where  ymeas  represents  the  “measured”  angle  6  and  a  represents  some  initial  angular 

quantity.  The  measured  values  are  then  perturbed  by  white  Gaussian  random  numbers  to 
simulating  sensor  noise  and  are  subsequently  fed  into  both  the  EKF  and  UKF.  Appendix 
A  -  Simple  Pendulum  Simulation,  shows  the  details  of  the  simulation,  including  the 
simulation  blocks,  and  associated  Matlab  code. 

F.  EKF  AND  UKF  ESTIMATION  RESULTS  USING  THE  SIMPLE 
PENDULUM  PROBLEM 

The  results  of  this  estimation  problem  show  how  the  EKF  does  not  estimate 
accurately  for  nonlinear  problems.  For  the  first  simulation  the  pendulum  was  set  to  0  = 
30°,  9=  0°/sec,  and  Iy  =  5  kg  m2.  Figure  9  shows  the  3a  plot  for  the  angle  error  between 
the  estimated  values  and  the  true  value.  The  3a  plot  is  typically  used  to  the  confidence 
interval  of  a  given  set  of  data.  While  the  term  “3a”  actually  refers  to  three  times  the 
variance  of  the  data  distribution,  mathematically  3a  can  be  translated  to  mean  that  our 
data  falls  within  approximately  99.73%  of  the  symmetric  confidence  interval  (Cl). 
Conversely,  this  means  that  approximately  0.27%  of  the  data  falls  outside  the  CL  The 
calculation  for  3  a  is  shown  below  where  the  variance  of  diagonal  values  of  the 
covariance  matrix  Pxx  are  used  for  n  number  of  states. 


Cl,l  Cl,2 

•••  cu 

C2,l  C2,2 

•••  C2,n 

2.59 

Cn,l  Cn,2 

Cn,n  \ 

:VC«(4)> 

i  =  1,2,... 

,,n 

2.60 

Here  we  can  see  that  the  EKF  cannot  accurately  estimate  the  state  due  to  the 
nonlinearity  of  the  system.  As  the  pendulum  swings  and  the  angle  increases,  the 
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nonlinearity  of  the  dynamics  increase  and  thus  the  filter  becomes  inaccurate.  Conversely, 
we  see  in  Figure  10  that  the  UKF  accurately  estimates  the  state  well  between  the  3a 
bounds. 
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Figure  9.  Angular  Errors  in  EKF  with  3a  Error  Bounds  Simulation  1  (large  0) 
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UKF  Angle  Error  (Deg) 


Figure  10.  Angular  Errors  in  UKF  with  3a  Error  Bounds  Simulation  1  (large  0) 


Figure  1 1  shows  similar  results  for  the  estimation  of  the  angular  velocity  co. 


EKF  Angular  Rate  Error  (Deg/Sec) 


Figure  11.  Angular  Rate  Errors  in  EKF  with  3a  Error  Bounds  Simulation  1  (large  0) 
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Figure  12.  Angular  Rate  Errors  in  UKF  with  3a  Error  Bounds  Simulation  1  (large  0) 

Furthermore,  we  can  see  that  the  estimation  for  the  moment  of  inertia,  Iy,  is  accurate  for 
both  EKF  and  UKF.  We  can  conclude  that  this  is  largely  because  Iy  is  a  constant 
quantity. 
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EKF  Moment  of  Inertia  Error 


Figure  13.  Moment  of  Inertia  Errors  in  EKF  with  3a  Error  Bounds  Simulation  1 

(large  0) 


UKF  Moment  of  Inertia  Error 


Figure  14.  Moment  of  Inertia  Errors  in  UKF  with  3a  Error  Bounds  Simulation  1 

(large  0) 
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These  plots  clearly  show  how  the  UKF  provides  a  more  accurate  solution  for  even 
simple  nonlinear  problems.  To  further  verify  this,  a  second  simulation  was  performed 

using  smaller  initial  conditions.  Using  6  =  \° ,  0  =  0,  and  I  =  5  kg  m2 ,  we  can  see  the 

both  filters  estimate  well  within  the  3a  bounds.  This  can  directly  be  associated  with  the 
small  angle  approximation  where  sin  0  =  0  for  sufficiently  small  angles.  These  plots  are 
shown  below. 


EKF  Angle  Error  (Deg) 


Figure  15.  Angular  Errors  in  EKF  with  3a  Error  Bounds  for  Simulation  2  (small  0) 
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UKF  Angle  Error  (Deg) 


Figure  16.  Angular  Errors  in  UKF  with  3a  Error  Bounds  for  Simulation  2  (small  0) 


Similarly,  as  shown  below,  we  can  see  that  the  angular  rates  also  fall  within  the 
bounds. 


EKF  Angular  Rate  Error  (Deg/Sec) 


Figure  17.  Angular  Rate  Errors  in  EKF  with  3a  Error  Bounds  for  Simulation  2 

(small  0) 
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UKF  Angular  Rate  Error  (Deg/Sec) 


Figure  18.  Angular  Rate  Errors  in  UKF  with  3a  Error  Bounds  for  Simulation  2 

(small  0) 


We  can  also  see  that  the  both  moments  of  inertia  also  converge  the  proper  values. 


EKF  Moment  of  Inertia  Error 
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Figure  19.  Moment  of  Inertia  Errors  in  EKF  with  3a  Error  Bounds  for  Simulation  2 

(small  0) 
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UKF  Moment  of  Inertia  Error 
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Figure  20.  Moment  of  Inertia  Errors  in  UKF  with  3a  Error  Bounds  for  Simulation  2 

(small  0) 
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III.  IMPLEMENTATION  OF  EKF  AND  UKF  FOR  SPACECRAFT 

ATTITUDE  DETERMINATION 

A.  GENERALIZATIONS 

While  Chapter  II  discusses  the  fundamentals  of  both  EKF  and  UKF,  this  chapter 
describes  the  implementation  of  both  methods  for  attitude  determination.  In  order  to 
perform  Kalman  filtering  for  attitude  estimation  we  must  first  examine  the  nature  of 
quaternion  estimation.  The  following  discusses  the  analytical  modeling  setup,  basic 
quaternion  attitude  kinematics,  and  finally,  the  implementation  of  both  EKF  and  UKF 
filters  for  spacecraft  attitude  estimation. 

B.  ANALYTICAL  MODELING  AND  SETUP  FOR  ATTITUDE 

DETERMINATION  SIMULATIONS 

1.  Background 

To  implement  the  Kalman  filters,  a  spacecraft  simulation  was  created  in 
MATLAB  Simulink.  Much  of  the  initial  foundation  for  this  simulation  was  built 
previously,  and  is  documented  in  [14]  and  [4],  For  a  better  understanding  of  how  the 
simulation  works,  the  following  sections  will  briefly  discuss  the  several  of  the  general 
Simulink  Blocks.  For  our  purposes,  we  will  define  the  general  simulation  blocks  as  the 
following 

•  Orbit  Propagator 

•  Environmental  Effects 

•  Dynamics  and  Kinematics 

•  Disturbance  Torques 

•  Sensors  and  Noise  Modeling 

Particularly,  the  blocks  that  will  be  discussed  are  Spacecraft  Kinematics,  Attitude 
Disturbance  Torques,  and  Spacecraft  Sensor/Noise  Modeling.  These  blocks,  and  the 
entire  simulation  are  shown  in  Appendix  B  -  Spacecraft  Attitude  Determination 
Simulation. 
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2.  Dynamics  and  Kinematics 

Shown  in  Figure  21,  the  Dynamics  block  calculates  the  spacecraft  angular  body 
rates  along  each  body  axis  by  integrating  applied  forces,  including  control  torques,  based 
on  Euler’s  equations  [15]. 


Figure  21.  Attitude  Dynamics  and  Kinematics  Simulink  Block 


The  Euler  equations  are  listed  below  [15]. 


co,  = 


Tx-(J7-Jv)(o7(ov 


CO , 


J, 

_  T  -(Jx-J_)ojxox 


CO,  =  - 


T-_-(Jy-Jx)®y®x 


J, 


2.61 


These  angular  rates  were  then  integrated  in  to  the  Kinematics  block  to  determine 
the  spacecraft  orientation.  For  these  simulations,  the  orientation  is  described  in 
quaternions.  The  quaternion  kinematic  differential  equation  is: 
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qA  0  0)z  -coy  cox  [qx 

q2  =  J_  -co,  0  ax  coy  q2 

q3  2  coy  ~cox  0  coz  q3 

q4  \  -My  0  Jk 

3.  Spacecraft  Attitude  Disturbance  Torques 

Three  major  torque  disturbances  were  taken  into  consideration  for  this  simulation, 
gravity  gradient,  aerodynamic  torque,  and  solar  torque.  For  gravity  gradient  torque,  the 
following  equations  were  used. 

(J,  ~Jy)C 2C3 
(  J  \  ~  *^y)ClC3 
y  ~Jx)C\C2 

where 


The  aerodynamic  torque  was  calculated  using  the  simple  drag  equation. 


2.64 


Faero=~pV2CDA 


Here  we  should  note  that  the  velocity  was  determined  from  the  orbit  parameters,  the 
coefficient  of  drag  (Cd)  was  assumed  to  be  2,  and  effective  area  (A)  was  determined 
using  spacecraft  component  areas  and  centers  of  pressure.  Finally,  solar  torque  was 
determined  using  a  Simulink  block  diagram  shown  in  Appendix  B. 
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4.  Spacecraft  Sensor  and  Noise  Modeling 

All  of  sensor  modeling  done  for  this  simulation  was  completed  in  Reference  [4], 
where  the  author  accurately  modeled  sensor  sampling  rates  and  noise  sources  based  on 
manufacturer  specifications.  We  can  see  this  in  the  “Sensors  Block”  of  the  simulation. 
The  most  important  information  from  the  previous  work  is  shown  in  Table  6,  which 
shows  the  example  of  noise  coefficients  for  the  gyro.  These  numbers  are  implemented 
into  the  spacecraft  simulation  gyro  random  noise  modeler. 


Table  6.  Summary  of  Gyro  Noise  Coefficients,  from  [4] 


£[?/;]  (°/Vsec) 

E\r,l ]  (°/Vsec3) 

Gyro  Data  1 

7.840e-04 

1.440e-07 

Gyro  Data  2 

7.840e-04 

3.240e-08 

Gyro  Data  3 

7.840e-04 

3.240e-08 

The  equation  for  modeling  the  internal  measurement  unit  (IMU)  is  listed  below. 

2.66 

2.67 


d)(t)  =  oit)  +  P(t)  +  riv(t) 


where  cb{t)  is  the  continuous-time  measured  angular  rate,  and  r/v(t)  and  r/u(t)  are 
independent  zero-mean  Gaussian  white-noise  processes. 
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C.  ATTITUDE  KINEMATICS  FOR  QUATERNION  ESTIMATION 

This  section  describes  the  Kalman  filter  as  it  applies  to  attitude  estimation.  It  is 
important  to  note  that  the  equations  found  in  this  section  apply  to  both  the  extended  and 
the  unscented  Kalman  filters. 


The  quaternion  is  defined  in  Equations  2.68,  2.69,  and  2.70. 

q  =  [sr 

s  =  [<h  <h  93f  =£sin(f) 


2.68 

2.69 


q4  =  cos(f)  2.70 

where  q  is  the  quaternion,  e  is  the  Euler’s  axis,  0  is  the  Euler’s  angle,  and  the  quaternion 
follows  the  normalization  of  qTq  =  1 .  The  attitude  matrix  can  be  related  to  the  quaternion 
by  the  equation  below  [2], 

^(q)  =  ST(q)'E(q)  2.71 


where  Sr(q)  and  xPr(q)  are  defined  by  Equations  2.72  and  2.73. 


Sr(q)  = 
^r(q)  = 


#4^3x3  +  [*?  X] 


^4^3x3  X] 


2.72 

2.73 


Here  /3x3  is  a  3x3 identity  matrix  and  [qx]  is  the  cross  product  matrix  described  below. 


41 


-<h 

0 


q  2 

-qi 

o 


2.74 


[«x] 


0 

qz 

_-q2 


qx 


Notably,  the  quaternion  error  cannot  accurately  found  by  subtraction,  as  the  result  would 
not  satisfy  the  unit  norm  constraint,  and  a  renormalization  would  be  needed.  The 
multiplicative  error  is  defined  as  [16]: 

S<\  =  q  <8>  q_1  2.75 

Here  we  use  the  symbol  <£>  to  indicate  the  quaternion  multiplication  [2],  This  relationship 
is  described  in  Equation  2.76. 


A(q'  )A(q)  =  A(q'  0q)  2.76 

For  implementation,  the  function  XI  was  used  in  Matlab.  This  can  be  seen  in  Appendix 
A  -  Matlab  Code  and  Simulink  Diagrams.  The  time  derivative  of  the  quaternion  error 
becomes 


£q  =  q®q1  +  q®q1  2.77 


As  derived  in  [16],  the  estimated  quaternion  kinematics  equation  is  given  by 


q(0=Ts[q(0]«)(0 


2.78 


Where  o>(t)  is  the  3x1  angular  velocity  vector.  The  local  error-quaternion, 
Sq=[sgT  Sq4J ,  can  now  be  used  to  find  the  generalized  Rodriguez  parameter  which  will 

be  useful  later  in  the  implementation  of  the  UKF  [4], 
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2.79 


8p  =  f 


*g 

a  +  Sq4 


where  a  is  a  parameter  from  0  to  1 ,  and  /  is  a  scale  factor.  Suggested  values  for,  /,  is 

2(a+l )  so  that  ||<Sp||  is  equal  to  &  for  small  errors,  where  3  is  the  angle  of  rotation  [4], 

In  the  simulations  presented  in  the  thesis  a  was  set  to  1  to  reproduce  the  results  shown  in 
[4],  While  the  propagation  of  the  state  and  covariance  can  be  accomplished  by  using 
numerical  integration  techniques,  the  measurement  observations  are  typically  sampled  at 
higher  rates  than  they  are  updated.  This  proves  useful,  as  we  can  use  a  discretized 
version  of  the  propagation  equations.  Using  the  power  series,  we  can  derive  the  new 
discretized  propagation  equation  from  2.78  [9]. 


e 


2 


2 


\2k 


Q  (co)t 


(2  *)! 


2k+\ 

(2k +  \) 

\ 

2.80 


Using  the  identities  described  in  Equations  2.81  and  2.82,  we  can  substitute  them  into 
Equation  2.80. 

\k  II  /v \\2k 


n2k(6?)=(-iy\\d)f  i 


4x4 


2.81 


n2*+1(®)  =  (-i)*||®||2*n(<») 


2.82 


Q.(co)t 


2  =/  Y 

1 4x4 

k= 0 


(-o' 


m\t 


2k 


(2*)! 


■+  a 


00 


(-1/ 


cot 


2k+\ 


k= 0 


(2^  +  1)! 


2.83 


This  equation  then  simplifies  using  trigonometric  identities  to  the  following 
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Cl(co)t 

e  2  =I4x4c os 


-\\co\\t 

V2y  j 


+ 


f  1 1.  A  I.  ^ 

sin  —  \\a>  t 

U"  11 J 

m\ 


2.84 


Finally,  the  quaternion  propagation  is  found  to  be  [9]: 


q-+i=Q(4+)q;  2.85 

Where  cb+k  and  i\k  are  the  post-update  estimates  and  12(7^  j  are  given  by  Equations  2.86 
and  2.87. 


2.86 


sin 


A 


cat  At 


0), 


J 


2.87 


For  both  the  EKF  and  the  UKF,  we  will  use  a  rate  gyro  and  a  magnetometer.  Given  a 
post-update  estimate  for  the  bias  fy ,  we  will  use  the  following  equation  to  find  the  post¬ 
update  angular  velocity  and  propagated  bias. 


^  +  ~  n+ 

°>k  =®k-Pk 

h*=k 


2.88 
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We  will  use  these  equations  in  both  the  EKF  and  UKF  to  solve  the  attitude  determination 
problem  in  the  following  sections. 

D.  CRASSIDIS  AND  MARKLEY’S  UNSCENTED  QUATERNION 

ESTIMATOR  (USQUE) 

In  this  section,  the  unscented  Kalman  filter  described  in  Crassidis  and  Markley’s 
paper  on  spacecraft  attitude  estimation,  Reference  [4],  is  implemented.  This  filter  is 
called  the  UnScented  QUatemion  Estimator  (USQUE).  More  specifically,  the  following 
describes  how  the  USQUE  is  implemented  in  spacecraft  attitude-determination 
simulations  using  MATLAB. 

First,  however,  we  must  take  step  back  and  look  at  implementation  of  the  UKF 
process  as  a  sequential  series  of  steps.  Figure  22  shows  the  UKF  graphically  in  the  form 
of  a  flow  chart.  Similarly,  we  will  refer  to  this  flow  chart  throughout  this  section  as  it 
follows  Crassidis  and  Markley’s  USQUE  closely.  The  Matlab  Code  generated  for  this 
simulation  also  follows  this  flowchart  and  is  listed  in  Appendix  A. 
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•  Calculate  Variance  for  Sensors 

•  Assume  initial  zero  bias 

•  Use  Sensor  information  for 
quaternion  and  rate  initialization 


•  Use  Cholesky  Decomp  to  Find  First 
Sigma  Points 

•  Find  the  Error  Quaternion  (Ref  [4] 
Eq.  33) 

•  Propagate  the  Updated  Quaternion 
(Ref  [4]  Eq.  34) 

•  Propagate  the  Updated  Bias  (Ref  [4] 
Eq.  35) 

•  Calculate  Mean  Points  (Ref  [4]  Eq.  7 
&  9,  Ref  [9]  Eq.  5.290  &  5.292) 


•  Calculate  Covariances  Terms  Pxx,  Pyy, 
and  Pxy  (Ref  [4]  eqs.  8,11,  and  1 3) 

•  Calculate  Gain  (Ref  [4]  4) 


•  Error  Covariance  Update 

•  State  Update 

•  Calculate  the  Updated  Error 
Quaternion  (Ref  [4]  45) 

•  Calculate  the  Final  Updated 
Quaternion  (Ref  [4]  44) 


•  Finally  the  error  quaternion,  del_p,  is 
set  to  zero  and  the  process  is  “reset” 
until  the  next  sensor  update. 


Figure  22.  Unscented  Quaternion  Estimator  Flow  Chart 
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To  begin,  Figure  23  shows  the  UKF  block  of  the  attitude  determination  simulation 
shown  in  Appendix  A. 

Table  7  describes  the  inputs  and  outputs  of  this  block. 


[w_  B  N  m  ]\> — ^ 

w_BN 

w_BNf_u 

- ►^W  BNf  u] 

q_BNf_u 

- *<1q  BNf  u] 

[Bm] 

Bm 

bias_f_u 

- ^<^[bias  f  u] 

Pdiag_u 

- ^IPdiag  u] 

~]  ► 

b 

Pnorm_u 

- ^<^fPnorm  u] 

Unscented  Kalman  Filter 


Figure  23.  Unscented  Kalman  Filter  Block  -  Level  1 


Table  7.  Description  of  Inputs  and  Outputs  for  UKF  Block  -  Level  1 


Input 

Output 

Variable  Name 

Description 

Variable  Name 

Description 

w_BNm,  ( cok ) 

Sensor  measured  angular 

rate  (from  gyro) 

w_BNf_u,  ( 6)  ) 

Estimated  angular 

rate 

Bm  ,(pk) 

Sensor  measured  magnetic 

field  (from  magnetometer) 

q_BNf_u,  (q ) 

Estimated  quaternion 

b,(A) 

Estimated  magnetic  field 

from  environment  model. 

bias_f_u 

Estimated 

magnetometer  bias 

Pdiagu 

Diagonal  terms  of  the 

covariance  matrix 

Pnormu 

Norm  of  the 

covariance  matrix 
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As  Figure  23  shows  the  top  level  of  the  UKF,  Figure  24  shows  the  Level  2  block 
showing  a  few  more  inputs.  These  will  be  discussed  further  in  the  section.  As  a  side 
note,  the  following  sections  are  also  well  documented  in  the  embedded  Matlab  code 
associated  with  this  Simulink  block.  This  code  can  be  found  in  Appendix  B. 


norm 


Figure  24.  Unscented  Kalman  Filter  Block  -  Level  2 


1.  Initialization 

Referring  to  Figure  22,  we  can  see  that  the  USQUE  process  begins  with  the 
initialization  portion  of  the  estimation.  Here  we  must  choose  the  initial  values  for  our 
states,  which  include  the  quaternion  and  the  bias.  For  the  initial  simulations,  the  initial 
quaternion  was  set  to  [0,0,0, 1]  and  the  bias  to  [0,0,0].  For  later  simulations,  as  discussed 
in  the  results  section,  initial  conditions  were  changed  to  highlight  major  differences 
between  the  UKF  and  EKF. 

2.  Calculation  of  Sigma  Points 

Calculations  of  Sigma  Points  begin  with  defining  the  following  state  vector: 
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zk(0)=*;  = 


k . 


2.89 


Here  we  can  use  Equation  2.79  for  8\>+k ,  which  is  the  4  x  1  error  quaternion,  and  the  3 

state  bias  term,  fy .  These  values  will  be  propagated  and  used  to  update  the  final  nominal 

state.  This  resulting  covariance  matrix  is  a  6  x  6.  It  is  important  to  note  here  that  for 
propagating  these  values  forward  we  can  now  use  Equations  2.42  through  2.50. 
However,  before  we  use  these  equations,  we  must  partition  the  sigma  points  %k(i)  so  that 
we  can  work  only  with  the  quaternion  portion. 


X*(  0  s 


Xtpd) 

_  xf  (0  _ 


/  =  0,1, _ ,12 


2.90 


where  yj/’  is  the  attitude  error  part,  and  /f  (i)  is  the  gyro  bias  part.  Now  that  we  have 

parsed  out  these  terms,  we  must  determine  the  new  quaternion  generated  by  multiplying 
the  error  quaternion  by  its  current  estimate. 

qi(0)  =  qi  2.91 


ql(0  =  <*q;(0®ql  *  =  f2,. ..,12 


2.92 


where  is  the  current  quaternion  estimate,  and  8<{+k  is  the  error  quaternion.  The  error 
quaternion  is  broken  up  into  the  3  state  quaternion  vector  8q+k ,  and  the  forth 
quaternion,  SqA/ ,  shown  in  Equations  2.93,  2.94,  and  2.95. 


Sq+k(i)  =  [Sq+kT(i)  dqlii) 


i  =  1,2, ..,12 


2.93 
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,  i  =  1,2,.. .,12 


2.94 


Sq+4t(i)  = 


-a  | 

|zr<of+N 

/2  +  (  1 

-a2) 

xf  (0  2 

f2  + 

X?(0 

2 

K(i)  =  r1  [>  +  <?<(0]x?(0,  »•  =  1,2,.  ..,12 


2.95 


We  chose  /  =  2  (a  + 1) ,  where  a  values  were  selected  using  Table  1  from  [4],  Next,  these 
updated  quaternions  are  propagated  forward  using  Equation  2.85  for  each  i,  or  step, 
shown  below. 

q-+1(0  =  Q(4+(0)q;(0  i  =  0,1,. ..,12  2.96 

where  again  the  angular  velocities  are  given  by  Equation  2.97  similar  to  Equation  2.88  in 
the  previous  section.  Here,  we  can  see  that  for  a>l  (0)  =  cbk  -/f  (0) ,  /f  (0)  =  PI . 

op(i)  =  cbk  (/),  i  =  0,1,. ..,12  2.97 

The  propagated  error  quaternions  are  now  calculated  using  Equation  2.98. 

^q*+,(0  =  q*+,(0  ®  [q*+,(0)]_1 ,  t  =  0,1 . 12  2.98 

it  is  interesting  to  note  that  where  c>q^+1(0)  here  should  be  the  identity  quaternion 

[0,  0,  0,  1],  Finally,  the  propagated  sigma  points  can  be  calculated  using  the  following 
equations. 


zZ(P)=o 


2.99 


xg,(o=/  <•= 1,2...., 12 

ci  +  Sq.Ji) 


2.100 


where  Sqk+i  (i)  and  Sq4  (z)  are  found  from  the  following  equation. 
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<?q*+i(0=[<%w(0  S(iZSl) 


2.101 


3.  Covariance  and  Gain  Calculations 

The  next  step  in  the  UKF  process  is  to  calculate  the  covariances  and  gains  which 
is  the  most  notable  difference  between  the  EKF  and  UKF.  Now  that  we  have  calculated 
our  sigma  points,  we  can  determine  these  values.  As  previously  mentioned  in  the 
derivation  of  the  UKF,  we  can  determine  the  predicted  covariance  matrix  P~+1 ,  shown  as 

Pxx  in  the  Matlab  code,  the  output  covariance  P^. \ ,  and  the  cross  correlation  covariance 
Pklx .  These  equations  are  found  as  2.39,  2.43,  and  2.45,  respectively.  These  equations 
are  utilized  in  the  “Covariance  and  Gain  Calculations”  section  of  the  embedded  Matlab 
code  for  the  UKF.  For  initial  conditions,  Pxx  is  the  set  to  Qk ,  where  the  variations  for  the 

sensors  are  used.  The  following  equation  is  used  for  Qk . 


q  _  6  ^  )^3x3  ^3x3 

2  L  03x3  <t„2/3x3_ 

The  mean  observation  is  also  needed  to  calculate  the  covariance  terms. 

x4[q"  (/)]?! 
x4[q_(0]r2 

^[q"(0]>v 


2.103 


2.104 
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With  these  covariance  matrices  calculated,  we  can  now  determine  the  Kalman  gain  K 
from  the  equation  in  Table  5.  Unscented  Kalman  Filter,  from  [9].  This  equation  is  also 
shown  below. 


Kk=p?ipn 


2.105 


4.  Update  Routine  for  States  and  Error  Covariance 

After  the  gains  are  calculated,  states  and  error  covariances  must  be  updated.  First, 
the  error  covariance  is  updated  using  the  following  equation. 

Pk+=P~-KkPrKTk  2.106 


The  state  update  is  found  using 


=xk+Kkvk 


2.107 


Finally  we  can  update  the  quaternions  using  the  following  set  of  equations  where 


*7f+l 


’  Pk  ■  i 


^+i=ql+i®<b+i(°) 


2.108 


2.109 


Sql  = 


SpLi  +fpf  +(1"a  )  fiPw 


f+  8p+ 


k+\ 


2.110 


=  .r'[a+s^t,]sPii 


2.111 
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These  equations  are  very  similar  to  those  used  earlier  to  find  the  initial  error  quaternion. 
The  final  step  here  is  to  update  the  bias  using  Equation  2.88.  For  further  clarification,  the 
Matlab  code  references  the  equations  used  with  respect  to  [4],  This  unscented  Kalman 
filter  was  built  to  be  compared  with  the  extended  Kalman  filter.  The  EKF  and  UKF 
models  for  attitude  are  based  on  the  model  presented  in  [9]  and  [4]  where  the  state  vector 
is  represented  as  the  error  in  the  quaternion  and  generalized  Rodriquez  parameter 
respectively. 

E.  IMPLEMENTATION  OF  THE  EXTENDED  KALMAN  FILTER 

The  EKF  implemented  in  this  simulation  uses  many  of  the  equations  used  in 
previous  sections.  Similar  to  the  previous  section,  Figure  25  shows  a  flow  chart  of  the 
EKF.  By  comparison,  we  can  see  very  clearly  that  the  major  difference  in  the  EKF  is  the 
calculation  of  the  sensitivity  matrix,  which  is  the  inherent  linearization  processes 
associated  with  this  filter.  Much  of  the  information  on  the  derivation  of  the  EKF  is 
discussed  in  [4]  and  Table  8  shows  a  summary  of  the  EKF  equations.  A  complete  listing 
of  both  the  Matlab  Simulation  block  diagram  and  embedded  Matlab  code  are  shown  in 
Appendix  B. 
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Calculate  Variance  for  Sensors 
Assume  Initial  Zero  Bias 
Use  Sensor  Information  for  Quaternion 
and  Rate  Initialization 


Calculate  the  Estimated  Quaternion 
(Ref  [9]  7.18) 


Calculate  the  Covariance  (Ref  [9]  7.57, 
7.58,  7.59,  and  7.60) 


Calculate  the  Sensitivity  Matrix  (Ref  [9] 
7.43) 

Calculate  Gain  (Ref[4]) 


Update  Covariance  (Ref  [9]  5.44) 
Update  Bias  (Ref  [9]  7.46) 
Calculate  the  Updated  Error 
Quaternion  (Ref  [9]  7.44) 
Calculate  the  Final  Updated 
Quaternion  (Ref  [9]  7.487) 


Finally  the  error  quaternion,  .5*5a,  is 
set  to  zero  and  the  process  is  “reset” 
until  the  next  sensor  update. 


Figure  25.  Extended  Kalman  Filter  Flow  Chart 
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Table  8. 


Summary  of  EKF  Equations,  from  [9] 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 
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IV.  COMPARISON  OF  SIMULATION  RESULTS  USING  EKF 

AND  UKF  FILTERING  METHODS 

A.  SIMULATION  CONDITIONS 

In  this  section,  several  performance  comparisons  between  the  USQUE  and  EKF 
are  made  through  simulations  using  the  previously  discussed  spacecraft  model  and  the 
designed  EKF  and  UKF  filters.  Using  a  500-km  circular  orbit  the  simulation  time  was  set 
at  4,000  seconds.  The  attitude  determination  hardware  in  these  simulations  consisted  of  a 
gyroscopic  rate  sensor  and  a  three-axis  magnetometer  (TAM).  The  magnetic  field 
reference  model  uses  a  magnetic  dipole  approximation  as  previously  discussed. 
Furthermore,  these  sensors  were  characterized  in  previous  work  [4],  In  the  first 
simulation,  the  initial  attitude  error  was  set  only  to  30°,  while  the  attitude  rate  error  was 
set  to  0°/sec  in  all  axes.  A  second  simulation  was  run  using  an  initial  attitude  error  of  30° 
and  an  attitude  rate  error  of  30°/sec  in  all  axes. 

B.  SIMULATION  1  RESULTS 

The  following  shows  the  results  of  Simulation  1.  Figure  26  shows  the  attitude 
error  of  the  quaternion  for  the  EKF  estimator  with  3cr  bounds.  We  can  see  that  the  EKF 
takes  approximately  4,000  seconds  before  the  error  is  bounded.  Conversely,  we  can  see 
that  the  attitude  error  of  the  UKF  is  bounded  in  approximately  2,500  seconds.  This  is 
shown  in  Figure  27  where  the  generalized  Rodriguez  parameters  are  shown.  It  is 
important  to  note  that  we  use  the  generalized  Rodriquez  parameters  instead  of  the 
quaternion  for  the  UKF  because  the  3cr  bounds  are  calculated  from  the  square  root  of  the 
diagonals  of  the  covariance  matrix  P .  For  the  UKF  the  covariance  matrix  is  built  from 
dp  which  is  shown  in  Equation  2.79,  as  dp  is  the  error  in  the  vector  of  the  generalized 
Rodriquez  parameter.  For  comparison  purposes,  Figure  30  will  show  both  normalized 
quaternion  errors  without  the  3cr  bounds. 
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Quaternion  Error  for  EKF 


Figure  26.  Simulation  1  Quaternion  Attitude  Error  with  3a  Bounds  for  EKF 
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Generalized  Rodriquez  Parameter  Error  for  UKF 


Figure  27.  Simulation  1  Generalized  Rodriguez  Parameter  Attitude  Error  with  for  3a 

Bounds  for  UKF 
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Similarly,  we  see  that  both  EKF  and  UKF  bias  errors  converge  in  a  similar  way. 
While  Figure  28  shows  the  bias  for  the  EKF  converging  within  the  3a  bounds  at 
approximately  2,700  seconds,  Figure  29  shows  convergence  at  a  little  after  1,000 
seconds.  We  can  also  see  that  the  initial  estimates  of  the  EKF  are  more  inaccurate  than 
the  UKF. 
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Figure  28.  Simulation  1  EKF  Bias  Errors  with  3a  Bounds 
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Figure  29.  Simulation  1  UKF  Bias  Errors  with  3a  Bounds 
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Finally,  we  can  see  that  the  normalized  EKF  and  UKF  attitude  errors  converge  as 
originally  predicted  and  demonstrated  in  the  simple  pendulum  problem.  Figure  30 
clearly  shows  the  better  performance  of  the  UKF.  Again,  as  a  nonlinear  estimator,  the 
UKF  consistently  shows  better  performance  on  all  figures.  This  is  again  shown  in 
comparison  plot  of  the  normalized  bias  errors  displayed  in  Figure  3 1  where  we  can  see 
that,  although  both  estimators  are  trending  appropriately,  that  the  UKF  performs 
significantly  better. 


Figure  30.  Comparison  of  EKF  and  UKF  Normalized  Attitude  Errors  for 

Simulation  1 
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Normalized  EKF  and  UKF  Bias  Errors 


Figure  31.  Comparison  of  EKF  and  UKF  Normalized  Bias  Errors  for  Simulation  1 

C.  SIMULATION  2  RESULTS 

The  second  simulation  shows  very  similar  results.  Although  we  can  see  similar 
trends  in  both  the  EKF  and  UKF  error  estimates,  we  can  see  that  the  UKF  consistently 
performs  better  in  every  plot.  Again,  Figure  32  shows  the  EKF  attitude  quaternion  error, 
which  settles  within  the  3a  bounds  at  approximately  3,000  seconds.  Figure  33  shows  the 
UKF  attitude  generalized  Rodriquez  parameter  error  settles  at  2,500  seconds.  The 
increased  performance  is  shown  without  fail  for  all  subsequent  plots  in  this  section. 
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Quaternion  Error  for  EKF 


Figure  32.  Simulation  2  Quaternion  Attitude  Error  with  3a  Bounds  for  EKF 
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Generalized  Rodriquez  Parameter  Error  for  UKF 


Figure  33.  Simulation  2  Generalized  Rodriguez  Parameter  Attitude  Error  with  for  3a 

Bounds  for  UKF 


Figures  34  and  35  show  that  the  bias  for  the  EKF  settles  at  approximately  2,700  seconds 
while  the  UKF  bias  settles  at  1,700  seconds. 
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Figure  34.  Simulation  2  EKF  Bias  Errors  with  3a  Bounds 
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Figure  35.  Simulation  2  UKF  Bias  Errors  with  3a  Bounds 
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Finally,  we  can  see  in  Figure  36  that  the  normalized  attitude  error  of  the  UKF  is 
much  better.  Similarly,  this  is  shown  in  Figure  37  with  the  comparison  of  the  normalized 
bias  errors. 


Figure  36.  Comparison  of  EKF  and  UKF  Normalized  Attitude  Errors  for 

Simulation  2 
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Normalized  EKF  and  UKF  Bias  Errors 


Figure  37.  Comparison  of  EKF  and  UKF  Normalized  Bias  Errors  for  Simulation  2 

D.  DISCUSSION  OF  RESULTS  AGAINST  PREVIOUS  LITERATURE 

Much  of  the  work  on  the  UKF  was  researched  [4],  In  this  paper,  Crassidis  and 
Markley  discuss  the  performance  of  the  UKF  as  it  applies  to  the  spacecraft  attitude 
determination  problem.  Figure  38  is  taken  from  [4]  and  shows  many  similarities  to  the 
spacecraft  model  designed  for  this  thesis.  Although  the  results  are  not  identical,  they 
shown  very  similar  trends  and  performance  characteristics.  It  should  be  noted  that  the 
initial  error  conditions  used  for  simulations  in  [4],  were  much  larger  and  are  used  here  to 
highlight  the  differences  in  accuracy.  The  simulations  done  for  this  thesis  were  set  at 
4,000  seconds. 
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V.  CONCLUSION 


A.  SUMMARY 

The  results  from  the  simulations  clearly  show  that  the  UKF  developed  here  is 
more  accurate  than  the  EKF  [2],  Both  the  EKF  and  UKF  were  rigorously  tested  and 
validated  against  previous  research  papers.  These  results  show  both  that  the  UKF  is 
largely  better  for  nonlinearities,  but  that  the  EKF  performs  rather  well.  To  take 
advantage  of  the  UKF,  large  nonlinearities  must  be  present  in  the  physical  dynamics  of 
the  system.  In  summary,  we  have  shown  that  the  UKF  has  a  lower  expected  error  than 
the  EKF  for  all  instances  of  spacecraft  attitude  determination.  We  showed,  in  the 
pendulum  problem,  that  as  the  nonlinearity  of  the  dynamics  increase,  that  the  UKF  shows 
increased  performance  over  the  EKF.  However,  for  slightly  nonlinear  or  linear 
estimation,  the  EKF  performs  well  and  will  provide  accurate  solutions.  The  one 
remaining  question  is  the  computational  expense  that  the  extra  computations  cost.  In  our 
simulations  the  UKF  performed  approximately  2.4  times  slower  than  the  EKF,  which  was 
consistent  with  [4],  As  the  optimization  of  any  process  is  measured  by  a  cost  function, 
one  must  evaluate  and  prioritize  the  resources  available.  Literature  tells  us  that  the  UKF 
has  2.5  times  the  cost  in  computational  time  of  the  EKF  [4],  For  spacecraft  with  relaxed 
attitude-control  requirements  and  low  computational  power,  it  could  be  argued  that  the 
EKF  could  perform  sufficiently  without  the  added  expense.  However,  the  UKF  can 
certainly  be  used  in  the  worst  case  conditions,  such  as  partial  loss  of  attitude  control,  and 
in  the  “lost  in  space”  scenario  where  anomalies  in  the  separation  event  from  the  launch 
vehicle  imparts  a  large  torque  on  the  spacecraft  hurling  it  into  an  unwanted  orientation. 
These  scenarios,  although  somewhat  unlikely,  mostly  likely  cannot  be  recovered  from 
using  an  EKF  estimator. 

B.  FUTURE  WORK 

This  thesis  presents  partial  validation  of  the  UKF  and  EKF  estimators.  Although 
the  results  are  favorable  and  largely  resemble  other  research  papers,  a  more  realistic 
simulation  would  require  hardware.  Further  developments  in  the  model  can  also  be 
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applied.  A  high  order  magnetic  field  model  could  be  implemented  if  the  computing 
resources  were  available.  Most  importantly,  Monte  Carlo  simulations  should  be  run  to 
show  the  full  performance  characteristics  of  both  filters.  All  of  this  work  is  completely 
possible  for  further  thesis  students  and  laboratory  research  such  as  the  currently  being 
performed  in  the  Nanosatellite  Advance  Concepts  Laboratory. 
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APPENDIX  A  -  SIMPLE  PENDULUM  SIMULATION 


SIMULINK  Block  Diagram  -  Simple  Pendulum  Simulation 


Extended  Kalman  Filer 
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Model  Initialization  Parameters 


time_step=0.05; 

sigmanoise=le-2; 

R=sigmanoiseA2; 

Q=diag([0.0000 1,0.0000 1,0.001]); 


fncA.m 


function  deriv  =  fncA(x) 

m=50  ; 

g=9.81; 

1=0 . 1; 

deriv=zeros  (3,1) / 
deriv ( 1 , 1 ) =x  ( 2 , 1 ) ; 

deriv (2, 1) =- (m*g*l*sin (x (1) ) ) / (x (3) +m*lA2) 


fncC.m 


%#eml 

function  y  =  fncC (x) 
Bo=0 .5; 

alf=30*pi/180 ; 
y=Bo*sin (x (1, 1) +alf ) ; 
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Unscented  Kalman  Filter  -  Embedded  Matlab  Block 


%#eml 

function  [x_kl ,  Pxx_kl ]  =  UKF (x_k, Pxx_k, Y_meas ,  ts , Q, R, kappa) 
%  This  block  supports  the  Embedded  MATLAB  subset. 

%  See  the  help  menu  for  details. 

Dx=size (x_k, 1) ; 

Dy=size (Y_meas, 1) ; 

NSig=2*Dx+l; 

sig_x= (chol ( (Dx+kappa) *Pxx_k) ) ' ; 

ooooooooooooooooooooo 

x_sig_k=x_k*ones ( 1 , NSig) +[ zeros ( Dx, 1 )  sig_x  -sig_x] ; 

RK= zeros (Dx, 4) ; 
x_sig_kl=zeros (Dx,NSig) ; 
y_sig_kl=zeros (Dy,NSig) ; 

for  i=l : NSig 

RK (  :  ,  1 ) =  fncA (x_sig_k (  :  ,  i )  )  ; 

RK ( : , 2) =fncA (x_sig_k ( : , i) +l/2*ts*RK ( : , 1) ) ; 

RK ( : , 3 ) =f ncA (x_sig_k ( : , i) +l/2*ts*RK (:,2)); 

RK ( : , 4 ) =  fncA (x_sig_k (  :  ,  i ) +ts*RK (  :  ,  3 )  )  ; 
x_sig_kl ( : , i) =x_sig_k ( : , i) +l/6*ts*RK* [ 1  2  2  1]'; 
y_sig_kl ( : , i) =  fncC (x_sig_kl ( : , i)  )  ; 

end 

W=ones (NSig, 1) / (2* (Dx+kappa) ) ; 

W  ( 1 , 1 ) =kappa/ ( Dx+kappa) ; 

x_klp=x_sig_kl*W; 

y_klp=y_sig_kl*W; 

Pxx_klp=Q; 

Pyy_klp=R; 

Pxy_klp=zeros (Dx,  Dy)  ; 
for  i=l : NSig 

xdif =x_sig_kl ( :  ,  i) -x_klp; 
ydif=y_sig_kl ( : , i) -y_klp; 

Pxx_klp=Pxx_klp+xdif *xdif 1 *W ( i , 1 ) ; 

Pyy_klp=Pyy_klp+ydif ^ydif ' *W ( i , 1 ) ; 

Pxy_klp=Pxy_klp+xdif *ydif ' *W ( i , 1 )  ; 

end 

K=Pxy_klp/Pyy_klp; 

Pxx_kl=Pxx_klp-K*Pxy_klp ' ; 
x_kl=x_klp+K* ( Y_meas-y_klp) ; 


75 


Extended  Kalman  Filter  Block 


Extended  Kalman  Filter  Block-Update  Block 


x  est  ukfl 
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Extended  Kalman  Filter  -Propagation  Block 


Update  of  X 
x_k+  K(h(x)-y)1 


77 


Integrate-  Embedded  Matlab  Block 


function  [xprop.  Phi]  =  Integrate (x, ts) 

%  This  block  supports  the  Embedded  MATLAB  subset. 
%  See  the  help  menu  for  details. 

KXl=fncA (x) ; 

KP1=F (x) ; 

KX2=f ncA (x+1/2 . 0*ts*KXl)  ; 

KP2=F (x+1/2 . 0*ts*KXl) * (eye (3) +1/2 . 0*ts*KPl) ; 
KX3=fncA (x+1/2 . 0*ts*KX2) ; 

KP3=F (x+1/2 . 0*ts*KX2) * (eye (3) +1/2 . 0*ts*KP2) ; 
KX4=fncA (x+ts*KX3)  ; 

KP4=F (x+ts*KX3) * (eye (3) +ts*KP3)  ; 

xprop  =  x  +  1/6 . 0*ts* (KX1+2*KX2+2*KX3+KX4) ; 

Phi  =  eye (3)  +  1/6 . 0*ts* (KP1+2*KP2+2*KP3+KP4) ; 
return 

function  deriv=F (x) 
m=50  ; 
g=9 .81; 

1= .  1  ; 

deriv=zeros  (3, 3) ; 
deriv ( 1 , 2 ) =1 ; 

deriv (2,1) =-m*g*l*cos (x (1, 1) ) / (x (3, 1) +m*lA2) ; 
deriv (2,3) =m*g*l*sin (x (1, 1) ) / (x (3, 1) +m*lA2) A2; 
return 


pendulum  plots.m 

x_true=zeros ; 

[m  n  p] =size  (x) ; 
x_true (l:m, 1 : p) =x ( 1 : m, l,l:p) ; 
x  true=x  true ' ; 


%%  State  1:  Theta 

%EKF  Error 
figure  ( 1 ) 

plot (t,  (x_est_ekf ( : , 1) -x_true ( : , 1) ) ^180/pi) ; 
grid  on 

title ( ' EKF  Angle  Error  (Deg) ') 

xlabel ( ' Time  (S) ' ) 

ylabel (’ Angle  Error  (Deg)') 

%  ylim ( [-5  5] ) 
hold  on 

plot (t,  3*sqrt (Pxx_est_ekf ( : , 1 ) ) * 180/pi ; 
plot (t,  -3*sqrt (Pxx_est_ekf (:, 1) ) *180/pi 1 , '-.'); 
%UKF  Error 
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figure ( 2 ) 

plot (t,  (x_est_ukf ( : , 1) -x_true ( : , 1) ) *180/pi) 
grid  on 

title ('UKF  Angle  Error  (Deg) ') 

xlabel ( 1  Time  (S)  '  ) 

ylabel (' Angle  Error  (Deg)') 

%  ylim ( [-5  5] ) 
hold  on 

plot ( t ,  3*sqrt (Pxx_est_ukf ( : , 1 ) ) * 180/pi ' , 
plot  (t,  -3*sqrt (Pxx_est_ukf ( : , 1) ) *180/pi ' ,  '  - 

%%  State  2 :  Omega  -  Angular  Rate 

%EKF  Error 
figure  ( 4 ) 

plot (t,  (x_est_ekf ( : ,  2 ) -x_true ( : , 2 ) ) *180/pi) 
grid  on 

title ( ' EKF  Angular  Rate  Error  (Deg/Sec) ') 
xlabel ( ' Time  (S) ' ) 

ylabel (' Angular  Rate  Error  (Deg/Sec) ') 

%  ylim ( [-10  10] ) 
hold  on 

plot  (t,  3*sqrt (Pxx_est_ekf ( :  ,  2 )  )  *180/pi  1  ,  '  -  . 
plot ( t ,  -3*sqrt (Pxx_est_ekf ( : , 2 ) ) *180 /pi ' ,  '  - 

%UKF  Error 
figure  ( 5 ) 

plot (t,  (x_est_ukf ( : , 2 ) -x_true (:,2))*180/pi) 
grid  on 

title ('UKF  Angular  Rate  Error  (Deg/Sec)') 
xlabel ( ' Time  (S) ' ) 

ylabel (' Angular  Rate  Error  (Deg/Sec) ') 

%  ylim ( [-10  10] ) 
hold  on 

plot  (t,  3*sqrt (Pxx_est_ukf ( :  ,  2 )  )  *180/pi  '  ,  '-. 
plot (t ,  -3*sqrt (Pxx_est_ukf ( : , 2 ) ) *180/pi ' , ' - 

%%  State  3:  Iy  -  Moment  of  Inertia 

%EKF  Error 
figure  ( 7 ) 

plot (t,  x_est_ekf ( : ,  3) -x_true ( : , 3) ) ; 
grid  on 

title ('EKF  Moment  of  Inertia  Error') 
xlabel ( ' Time  (S) ' ) 

ylabel (' Moment  of  Inertia  (kg  mA2) ') 

%  ylim ( [-5  5] ) 
hold  on 

plot (t,  3*sqrt (Pxx_est_ekf (:, 3) )','-.') ; 
plot (t,  -3*sqrt (Pxx_est_ekf ( : , 3 ))','-.') ; 

%UKF  Error 
figure  ( 8 ) 

plot (t,  x_est_ukf ( : ,  3) -x_true ( : , 3) ) ; 
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grid  on 

title ( ’ UKF  Moment  of  Inertia  Error') 
xlabel ( ' Time  (S) ' ) 

ylabel (' Moment  of  Inertia  (kg  m\A2) ') 

%  ylim ( [-5  5] ) 
hold  on 

plot (t,  3*sqrt (Pxx_est_ukf (:, 3 ) 
plot (t,  -3*sqrt (Pxx_est_ukf ( : , 3) ' 
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APPENDIX  B  -  SPACECRAFT  ATTITUDE 
DETERMINATION  SIMULATION 


ADS  SpacecraftSim.mdl 


— ►<^[J_Matrix] 
Interia  Matrix 


- ►  [w  BNo] 

Initial  Inertial  Body  Rates 


- -►*  [w_ON] 

Orbital  in  Intertial  Rates 


-KjmulJ 

Earth  Gravity  Const 


[bShat] 

[eclipse] 


[rho_a] 


[w_BNm] 


No  Init  Attitude  Err 
[0,0, 0,1]’ 


C_BO 

radius 

J_Matrix 

Tangential  Velocity 
Radial  Velocity 
-►  solar  vector 
-►  Eclipse 
-►  rho_a 


Disturbance  Torques 


Unscented  Kalman  Filter 
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Orbit  Propagator-Simulink  Block 


Calculate  Orbital 
Elements 


Calculate  Altitude-Simulink  Block 


Re 


Keplerian  Orbit  Propagation-  Simulink  Block 


R_eci 
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Calculate  Orbital  Elements  -  Embedded  Matlab  Code 


function  [ rho_earth, beta, TA, r , Vt , Vr ,  Lat ,  Long]  = 
fen (time, Re, incl,ecc,R,V) 

% 

o 

o 

%  This  function  computes  the  classical  orbital  elements 
%  from  the  state  vector  (R,V)  using  Algorithm  4.1.  As  well  as 
%  other  orbital  parameters  needed  by  the  model. 

o 

o 

%  mu  -  gravitational  parameter  (itT3/s~2) 

%  R  -  position  vector  in  the  geocentric  equatorial  frame  (m) 

%  V  -  velocity  vector  in  the  geocentric  equatorial  frame  (m/s) 
%  r,  v  -  the  magnitudes  of  R  and  V 
%  vr  -  radial  velocity  component  (m/s) 

%  H  -  the  angular  momentum  vector  (irT2/s) 

%  h  -  the  magnitude  of  H  (m~2/s) 

%  incl  -  inclination  of  the  orbit  (rad) 

%  N  -  the  node  line  vector  (m~2/s) 

%  n  -  the  magnitude  of  N 
%  cp  -  cross  product  of  N  and  R 

%  RA  -  right  ascension  of  the  ascending  node  (rad)  not  used 
************** 

%  E  -  eccentricity  vector 
%  ecc  -  eccentricity  (magnitude  of  E) 

%  eps  -  a  small  number  below  which  the  eccentricity  is 
%  considered  to  be  zero 

%  w  -  argument  of  perigee  (rad)  not  used 
******************************* 

%  TA  -  true  anomaly  (rad) 

%  Vt  -  tangential  velocity  (m/s) 

%  Vr  -  radial  velocity  (m/s) 

%  rho_earth  -  earth  anglular  radius 
%  beta  -  beat  angle  (rad) 

%  Lat  -  Latitude  of  satellite  (rad) 

%  Long  -  Longitude  of  satellite  (rad) 

o 

o 

mu  =  398 . 6004418el2;  %  mA3/sA2 

eps  =  1.0e-10; 
r  =  norm (R) ; 
v  =  norm (V) ; 
vr  =  dot(R,V)/r; 

H  =  cross (R, V) ; 
h  =  norm (H) ; 

%  Calc  inclination 

%{ 

c  =  H  (3)  /h; 

if  (c  <  -1)  &&  (c  >  1) 

c  =  c  -  pi; 

end 

incl  =  acos (c) ; 

%} 

%  Calc  right  ascension  of  the  ascending  node  (rad) 

N  =  cross ([0  0  1 ] , H ) ; 
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n  =  norm (N) ; 


%  Calc  Eccentricity 

E  =  1/mu* ( (vA 2  -  mu/r) *R  -  r*vr*V) ; 

%ecc  =  norm (E) ; 

%  True  Annomoly 
if  ecc  >  eps 

c  =  dot (E, R) /ecc/r  ; 
if  (c  <  -1)  &&  (c  >  1) 

c  =  c  -  pi; 

end 

TA  =  acos (c) ; 
if  vr  <  0 

TA  =  2*pi  -  TA; 

end 

else 

cp  =  cross (N, R) ; 
c  =  dot (N, R) /n/r ; 
if  (c  <  -1)  &&  (c  >  1) 

c  =  c  -  pi; 

end 

if  cp ( 3 )  >=  0 

TA  =  acos (c) ; 

else 

TA  =  2*pi  -  acos(c); 

end 

end 

%  Calculate  the  tangential  and  radial  velocities 
Vt  =  h/r; 

Vr  =  mu/h*ecc*sin (TA) ; 

%  Calculate  earth  angular  radius 
rho_earth  =  asin (Re/r) ; 

%  Beta  calcs 
wb_0  =  0 ; 
ub_0  =  0; 

wb_dot  =  (-9.9639/86400) *rho_earthA (3.5) *cos (incl) / (l-eccA2) A2; 

wb  =  (wb_0  +  wb_dot*time) *pi/180 ; 

gamma  =  23 . 442*pi/180 ;  %rad 

ub_dot  =  (0 . 985648/86400) *pi/180;  %rad 

ub  =  ub_0+ub_dot*time; 

beta  =  asin ( sin (ub) *sin (gamma) *cos ( incl )  +  ... 

cos  (ub) *sin (incl) *sin (wb) -sin (ub) *cos (gamma) *sin (incl) *cos (wb) ) 

%  For  Mag  Calc  - 

%  Calculate  Earth  Coordinate  by  Simulate  the  Earth's  Rotation 
%  Track  the  movement  of  (0  Lat,  0  Long) 

PhiE  =  time*2*pi/ (23. 93*3600) ; 

ThetaE  =  0; 
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EarthCoord  =  [Re, ThetaE, PhiE] ; 

%  Calculate  Sat  Coordinate  in  Polar 
X  =  R  ( 1 )  ; 

Y  =  R  ( 2  )  ; 

Z  =  R  ( 3 )  ; 

Theta  =  atan2 ( sqrt (XA2+YA2 ) , Z ) ; 

Phi  =  atan2 (Y,X) ; 

SatCoord  =  [ r, Theta, Phi ] ; 

%  Calculate  Lat  and  Long 
ThetaO  =  pi/2-SatCoord ( 2 ) ; 

PhiO  =  SatCoord(3); 

Phil  =  EarthCoord ( 3 ) ; 

Lat  =  ThetaO; 

Long  =  (PhiO-Phil) ; 

if  Long  >  pi 

Long  =  Long-2*pi; 

end 

if  Long  <  -pi 

Long  =  Long+2*pi; 

end 
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Environmental  Effects-  Simulink  Block 


— ►(  i  ) 

rho_a 

—►CO 

Magnetic  Field 

—►CO 

Body  Sun  Vector 

— ►CO 

Inertial  Sun  Vector 

-KJ> 

Eclipse 
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EFI  2  ECEF-Embedded  Matlab  Code 


function  [C_EN,C_NE]=  ECI_2_ECEF (EarthRotation,  time) 
%  transformation  of  eci  to  ecef  coordinates 


theta  =  EarthRotation ( 3 ) *time ; 


C_EN  =  [  cos (theta)  -sin  (theta)  0; 

sin  (theta)  cos (theta)  0; 


0 

C  NE  =C  EN ' ; 


0 


l] ; 


Magnetic  Dipole  Model  -  Embedded  Matlab  Code 

function  B_eci  =  Earth_Mag_Field (R) 

%  Magnetic  dipole  model  -  in  Tesla 


theta  =  11.7; 


o. 


deg 


DCM  =[10  0;  0  cosd(theta)  sind(theta);  0  -sind(theta)  cosd (theta) ] ' 


%  N/AmpA2 
%  A*mA2 


muO  =  4*pi*10e-7; 

M  =  DCM* [0  0  8e22 ]  '  ; 


r  =  norm (R) ; 
r_hat  =  R/r; 

B_eci  =  muO* (3*dot (M, r_hat) *r_hat-M) / (4*pi*rA3) ; 


S  hatl  -  Embedded  Matlab  Code 

function  [ S_body, S_inertial ]  =  S (C_BO, C_NB, beta, TA) 

%  This  block  supports  an  embeddable  subset  of  the  MATLAB  language. 
%  See  the  help  menu  for  details. 

B=beta; 

S_orbit  =  [cos (B) *sin (TA) ; . . . 

sin  (B)  ;  .  .  . 
cos (B) *cos (TA) ] ; 

S_body  =  C_BO*S_orbit ; 

S_inertial  =  C_NB*S_body; 
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Dynamics  and  Kinematics  -  Simulink  Block 


i 


disturbance  torque 

:  -f  j|  > - 

control  torque 


5 


Moments  of  Inertia 


Initial  inertial  body  rates 


omega  of  orbit  in  inertial 


Torques 

w_BN 

Moments  of  Inertia 

initial  inertial  body  rates 

w_ON 

omega  of  orbit  in  inertial 

C_BO 

omega_BN 

q_BN 

q_BO 

EulerB0321 

omega_ON 

EulerBN321 

C_NB 

Attitude  Dynamics 


Kinematics 


-xl 5  J 

w  BN 


#(1) 

C_BO 

►  2 

q_BN 


►  3  ! 

q_BO 


EulerBO  321 


EulerBN  321 

C  NB 


Attitude  Dynamics  -  Simulink  Block 


(  4  — 

omega  of 
orbit  in  inertial 


►GD 

w_ON 
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torque2omegadot.m  EML 


function  Wdot  =  torque2omegadot (T,  J,  W) 

%  This  function  takes  input  of  applied  torque  (T)  in  component 
%  elements,  current  angular  velocity  (W)  in  component 
%  elements,  and  the  moment-of-inertia  matrix  (J)  as  a  diagonal 
%  matrix  containing  the  MOIs  for  the  principal  axes  of  the  body 
%  along  the  diagonal.  Angular  acceleration  is  then  computed  and 
%  output  as  a  3x1  vector  (Wdot) . 

W x  =  W  ( 1 )  ;  Wy  =  W(2);  Wz  =  W(3); 

Jxx  =  J ( 1 , 1 ) ;  Jyy  =  J(2,2);  Jzz  =  J(3,3); 

Tx  =  T ( 1 )  ;  Ty  =  T ( 2 )  ;  Tz  =  T (3)  / 

wdotx  =  (Tx- ( Jzz- Jyy) *Wz*Wy) / Jxx; 
wdoty  =  (Ty- ( Jxx- Jzz ) *Wx*Wz ) / Jyy; 
wdotz  =  (Tz- ( Jyy- Jxx) *Wy*Wx) / Jzz ; 

Wdot  =  [wdotx;  wdoty;  wdotz]; 
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Disturbance  Torques-  Simulink  Block 


Torque 

Solar 
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Torque  Gravity  Gradient  -  EML 


function  Tgg  =  T_Grav_Grad (C_BO,  r,  J,  mu) 

%  T_Grav_Grad  takes  inputs  of  the  spacecraft  inertia  matrix  ( J) , 

%  current  orbit  radius  (r)  in  m,  and  the  Orbit-to-Body  Frame  DCM  (C_BO) 
%  to  calculate  the  gravity  gradient  torque  in  the  body  frame  (Tgg)  and 
%  orbit  frame  (T_o) .  The  orb_vec  vector  defines  which  orbit  frame  axis 
%  is  aligned  with  the  force  producing  the  torque.  In  this  case,  the  z- 
axis 

%  points  along  nadir  in  the  orbit  frame,  and  corresponds  to  the  r- 
vector 
%  direction. 

orb_vec  =  [0;  0;  1] ; 

c  =  C_BO*orb_vec; 

Jxx  =  J ( 1 , 1 ) ; 

Jyy  =  J (2, 2)  ; 

Jz  z  =  J  ( 3 ,  3 )  ; 

Tgg  =  3*mu/rA3* [ ( Jzz-Jyy) *c (2) *c (3)  ;  .  .  . 

(Jxx-Jzz) *c(l)*c(3);... 

( Jyy-Jxx) *c (1) *c (2) ] ; 


Torque  Aerodynamic  -  Simulink  Block 


Vr  frame  and  V  A2  Calucate  aero  torque 

by  element  side  and 
sum  for  total 
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Determine  Normalized  Velocity  Vector  -  Simulink  Block 


1 

C_BO 


2 

V  orbit 


magnitude 


Calculate  Aero  Torque  -  Simulink  Block 


Sum  of 
Elements 
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Solar  Torque  -  Simulink  Block 


1 

solar  vector 


Calculate  Solar 
Torque  from  side  elements 

2 - 

Eclipse 


0 


Switch 

m 


K  1 
Tsolar 


Calculate  Solar  Torque  -  Simulink  Block 


1 

s_body 


Sum  of 
Elements 


<  1  ) 

Tsolar 
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Sensors  -  Simulink  Block 
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Gyro  -  Simulink  Block 


Star  Tracker  -  Block  Diagram 


q_BNm 
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Star  Tracker  -  Embedded  Matlab  Code 


function  q  =  StarTracker (u , qbn) 

ph=u ( 1 ) / 2 ;  th=u(2)/2;  ps=u(3)/2; 

sph  =  sin  (ph) ;  sth  =  sin(th);  sps  =  sin(ps); 
cph  =  cos (ph) ;  cth  =  cos (th) ;  cps  =  cos  (ps) ; 

q  =  [ sph*cth*cps-cph*sth*sps ; 
cph*sth*cps+sph*cth*sps ; 
cph*cth*sps-sph*sth*cps ; 
cph*cth*cps+sph*sth*sps ] ; 

[Y  I ] =max (abs (q) ) ; 

q=q/ norm (q) *sign (q  ( I , 1 ) /qbn (1,1)); 


Sun  Sensors  -  Block  Diagram 

DCM  for 
Sun  Sensor  1 


ATT.m 

function  att  =  ATT (  quat  ) 

att  =  transpose (XI (quat) )  *  PSI (quat) ; 

return 


96 


Sun  Sensor  Facing  YO  -  Simulink  Block 


Sun  Sensor  Facing  YOl  -  Block  Diagram 

Ss 


Lookup 
Table  (2-D) 
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Mag  Inertial  to  Body  -  Simulink  Block 


MAGNETIC  FIELD  TRANSFORMATION  from  INERTIAL  to  BODY  COORDINATES 

Product  ATT 


B 

-ntx- 


Matrix 

Multiply 


◄ 

◄ 


att  ATT  quat 


q_BN 

T" 

b 

2 


Magnetometer  Model  -  Simulink  Block 


98 


Multiplicative  Quaternion  Extended  Kalman  Filter-  Simulink  Block 
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Bm 

' 3  D 

b 

t  ekf 


sig 


ST_on 

SS- 


1881  - 
rss2  - 


Gain  ^  q 


wklt 

wkl 

q_init 

bias_init 

qkl 

Bkl  EKF 

B 

biaskl 

dt 

Pkl 

sig 

Gain  1 
mag  _on 


ST_on 

SS1_on 

SS2_on 

Flag 

MGon 

Generate  Measurment 
Flag  Vector 


I  LA  It 
inctior 
diag 


5 1  LAt 
inctior 
norm 


-►  [Flag] 


w  BNf 


X  2 

q_BNf 


-X  3  ) 

bias  f 


>  4 

Pdiag 


►  5  ) 

Pnorm 


99 


MEKF  -  Embedded  Matlab  Code 


function  [wkl , qkl , biaskl ,  Pkl  ]  =  EKF (wklt ,  q_init , bias_init , Bkl , B, dt ,  sig) 

sig_v  =  sig ( 1 ) ; 
sig_u  =  sig (2) ; 
sig_mag  =  sig (5); 

persistent  qk  biask  wk  Pk; 

%  Initialize  States  and  Measurement 
if  isempty(qk) 
qk=q_init ; 
biask  =  bias_init; 
wk  =  wklt; 

Pk=[  (0.8) A2  *eye ( 3 )  zeros(3);  zeros (3)  (3*pi/180) A2*eye (3) ] ; 


wkl=wk; 

qkl=qk; 

biaskl=biask; 

Pkl=Pk; 

return; 

end 

%%  Propagation 
biaskl  =  biask; 


Skew_w  =  SKEW (wk) ; 

Mag_w  =  norm(wk) ; 

psik  =  (sin ( l/2*Mag_w*dt ) /Mag_w) *wk; 

Omega  =  [cos ( l/2*Mag_w*dt ) *eye (3) -SKEW (psik)  psik; 

-psik'  cos ( l/2*Mag_w*dt )  ]; 

qkl  =  Omega*qk; 

Phi_ll  =  eye ( 3 ) -Skew_w*sin (Mag_w*dt ) /Mag_w  +  Skew_wA2* fl¬ 
ees  (Mag_w*dt)  )  /Mag_wA  2  ;  %  7.59b 

Phi_12  =  Skew_w* ( 1-cos (Mag_w*dt) ) /Mag_wA2  -  eye(3)*dt  -... 

%  7.59c 

Skew_wA2* (Mag_w*dt-sin (Mag_w*dt) ) /Mag_wA3; 

Phi_21  =  zeros (3); 

%  7 . 5  9d 

Phi_22  =  eye (3); 

%  7 . 5  9e 

Phi  =  [ Phi_l 1  Phi_12 ;  Phi_21  Phi_22] ; 

%  7.59a 

Gk  =  [-eye (3)  zeros (3) ;  zeros (3)  eye (3)]; 

Qk  =  [  (sig_vA2*dt+l/3*sig_uA2*dtA3) *eye (3)  - ( l/2*sig_uA2*dt A2 ) *eye (3) 
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- (l/2*sig_uA2*dtA2) *eye (3) 


(sig_uA2*dt) *eye (3) 


]  ; 

Pkl  =  Phi*Pk*Phi ' +Gk*Qk*Gk '  ; 

%%  Update  Loop  - 

Att  =  ATT (qkl ) ; 
delX  =  zeros  (6,1); 

%  Update  for  Magnetometer  Measurement 
H  =  [ SKEW (Att*B)  zeros ( 3 , 3 )  ]  ; 

R  =  sig_magA2*eye (3)  ; 

%  Gain 

K  =  (Pkl*H' ) / (H*Pkl*H'  +  R) ; 

%  Update 

Pkl  =  (eye (6)  -  K*H)*Pkl; 

res  =  Bkl  -  Att*B; 

delX  =  delX  +  K* (res-H*delX) ; 

qkl  =  qkl+l/2*XI (qkl) *delX(l:3, :) ; 
qkl  =  qnormalize (qkl 1 *qkl r qkl ) ; 

biaskl  =  biaskl  +  delX(4:6,:); 

wkl  =  wklt  -  biaskl; 

%  Save  previous  values 
qk  =  qkl ; 
biask  =  biaskl; 
wk  =  wkl; 

Pk  =  Pkl; 

return 

o 

o 


%%  Normalizing  routine  for  quaternions 
function  qkl  =  qnormalize (qnorm, qkl ) 
while  (qnorm)  >  1 

if  qnorm  <  1  +  le-9 

qkl  =  ( (3  +  qnorm) /(I  +  3*qnorm) ) *qkl ; 
%  rescale  quaternion  to  (errA3)/32 

else 

qkl  =  qkl/sqrt (qnorm) ; 

%  renormalize  quaternion 

end 

qnorm  =  qkl!*qkl; 

end 

return 

o 

o 
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Unscented  Kalman  Filter  -  USQUE-  Simulink  Block 


►CO 

w  BNf  u 


q_BNf_u 


bias  f  u 


MATLAB 

Function 


diag 


MATLAB 

Function 

norm 


4 

Pdiag  _u 


5 

Pnorm  u 
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UKF-  USQUE-  Embedded  Matlab  Code 


function  [wkl , qkl , biaskl ,  Pxx_kl ]  = 

UKF (wklt ,  q_init , bias_init , Bkl ,  B , dt ,  sig, lambda, a) 

o 

o 


%  Initialization 

o 

o 


%  Variance  of  Sensors 
sig_v  =  sig ( 1 ) ; 

sig_u  =  sig (2 ) ; 

sig_mag  =  sig (5); 

f  =  2*  (a+1 ) /  %  Ref  [3]  pg  6 


persistent  qk  biask  Pxx_k 
%  Initialize  States  and  Measurement 
if  isempty(qk) 
qk=q_init ; 
biask  =  bias_init; 
disp (q_init ) ; 
disp (bias_init) ; 

Pxx_k= [  ( 1 ) A2 ^eye ( 3 )  zeros (3);  zeros  (3)  (3*pi/180) A2^eye (3) ] ; 

end 

o 

o 


o 

o 

o 

o 


Calculation  of  Sigma  Points 


Qbar_k  =  dt/2^ [ (sig_vA2-l/6*sig_uA2*dtA2 ) ^eye (3)  zeros (3)  ; 

%  Ref  [3]  42 

zeros (3)  sig_uA2 ^eye ( 3 )  ]; 


%  Sigma  points  equations 
x_k  =  [ [0  0  0] ’ ;  biask] ; 

Dx=size (x_k, 1) ; 

Dy=size (Bkl,  1)  ; 

NSig=2*Dx+l; 

sig_x=chol ( (Dx+lambda) * (Pxx_k+Qbar_k) ) 1 ;  %  Ref 

[3]_  5a 
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Ref  [3] 


chi_sig_k=x_k*ones ( 1 , NSig) +[ zeros (Dx, 1 )  sig_x  -sig_x] ; 

5  a 

chi_sig_kl=zeros ( 6 , NSig) ;  %  Ref  [3]  32 

y_sig_kl=zeros (3, NSig) ; 
q_kl=zeros (4,1) ; 


for  i=l : NSig 

del_q_k=delp2delq  ( chi_sig_k  ( 1 :  3 ,  i )  ,  a,  f )  ; 
q_sig_k  =  quaterr (del_q_k,  qk) ; 

o 

o 


%  Propagate  Forward  the  Quaternion  (still  in  the  loop!) 

o 

o 


w_sig_k  =  wklt  -  chi_sig_k ( 4 : 6 , i ) ; 

Mag_w  =  norm (w_sig_k) ; 

psik  =  (sin ( l/2*Mag_w*dt ) /Mag_w) *w_sig_k; 
zk  =  cos (l/2*Mag_w*dt) ^eye (3) -SKEW (psik) ; 

Omega  =  [  zk,  psik; 

-psik',  cos ( l/2^Mag_w^dt )  ]; 

q_sig_kl  =  Omega *q_sig_k; 

o  _ 

O - 

%  Saving  the  q ( — ) k+ 1(0) 

o  _ 

0 - 

if  i==l 

q_kl=q_sig_kl ; 

end 


%  Ref  [3]  35 


%  Ref  [3]  29 

%  Ref  [3]  34 

%  Ref  [3]  36 


del_q_kl  =  quaterr (q_sig_kl ,  [ -q_kl ( 1 : 3, 1 ) ; q_kl ( 4 , 1 ) ] ) ; 

chi_sig_kl (1 : 3, i)  =  f ^del_q_kl ( 1:3 , 1 ) / (a+del_q_kl ( 4 , 1 ) ) ;  %Ref  [3] 

37b 

chi_sig_kl ( 4 : 6 , i )  =  chi_sig_k(4:6,i); 
y_sig_kl ( : , i)  =  ATT (q_sig_kl) ^B; 

o 

o 


%  Note:  The  bias  does  not  change  so  chi_sig_kl ( 4 : 6, i)  stays  the 

same 

o 

o 


end 
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Following  USQUE  Method  Ref  [3]  pg  6 


%  Calculating  Weights 
R  =  sig_magA2*eye (3) ; 

W=ones (NSig, 1) / (2* (Dx+lambda) ) ; 
W ( 1 , 1 ) =lambda/ (Dx+lambda) ; 

%  Mean  Point  Calculations 

x_klp=chi_sig_kl*W; 

y_klp=y_sig_kl*W; 


%  Covariance  and  Gain  Calculations 

o 

o 


%  Error  Covariance  Calculation 
Pxx_klp=Qbar_k; 

Pyy_klp=R; 

Pxy_klp=zeros (Dx,  Dy)  ; 
for  i=l : NSig 

xdif =chi_sig_kl (  :  ,  i) -x_klp; 
ydif=y_sig_kl  (pi)  -y_klp; 

Pxx_klp=Pxx_klp+xdif *xdif ' *W ( i , 1 )  ; 
Pyy_klp=Pyy_klp+ydif *ydif ' *W ( i , 1 )  ; 
Pxy_klp=Pxy_klp+xdif *ydif 1 *W ( i f  1 )  ; 

end 

%  Gain  and  Update 
K  =  Pxy_klp/Pyy_klp;  %  Gain 

Pxx_kl  =  Pxx_klp-K*Pxy_klp ' ;  %  Error  Covariance  Update 

x_kl  =  x_klp+K* (Bkl-y_klp) ;  %  State  Update 

%  Calculation  of  Updated  Quaternion! 
del_q_kl=delp2delq (x_kl (1:3, : ) , a, f ) ; 
qkl  =  quaterr (del_q_kl , q_kl ) ; 
qkl  =  qnormalize (qkl )  ; 
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biaskl  =  x_kl (4 : 6,1) ; 
wkl  =  wklt-biaskl; 


%  Setup  for  Next  Update 

o 

o 


qk=qkl ; 
biask=biaskl ; 
Pxx_k  =  Pxx_kl; 

return 

o 

o 

o 

o 

o  o 
o  o 


%%  Normalizing  routine  for  quaternions 
function  qkl  =  qnormalize (qkl ) 
qnorm=qkl ' *qkl  ; 
while  (qnorm)  >  1 

if  qnorm  <  1  +  le-9 

qkl  =  ( (3  +  qnorm) /(I  +  3*qnorm) ) *qkl / 
%  rescale  quaternion  to  (errA3) /32 

else 

qkl  =  qkl/sqrt (qnorm) ; 

%  renormalize  quaternion 

end 

qnorm  =  qkl'*qkl; 

end 

return 

o 

o 
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o\o  o\°  o\o 


ADS  MainScript.m 


%%  Spacecraft  Attitude  Determination  Script 
%  Note  that  this  code  runs  both  the  EKF  and  the  UKF 


Created  by  Orlando  X.  Diaz 
Advisor  Dr.  Marcelo  Romano 
Co-Advisor  Dr.  Hyun-wook  Woo 

%%  Format 

clear  all 
close  all 
clc 

global  CONST 
R2D  =  180/pi; 

D2R  =  pi/180; 

%%  Set  Simulation  Conditions 

InitialEuler  =  [0,0,0] ;%deg 
Ref erenceEuler  =  [00  0];%deg 


%***  Toggle  switches  turn  the  labeled  functions  on  (1)  or  off  (0) . 


k  k  k 

Tgg_toggle  =  1 ;  % 

Taero_toggle  =  1;% 

Tsolar_toggle  =  1;% 

timeOn  =  1; 

taOn  =  0; 

cboOn  =  0 ; 

qbnOn  =  1 ; 

qbnmOn  =  1  ; 

rOn  =  0 ; 

hOn  =  0 ; 

e3210n  =  1; 

wbnOn  =  1 ; 

tcOn  =  0; 

hsOn  =  0; 

wbnfOn  =  1; 

biasOn  =  1; 

biasfOn  =  1; 

pdOn  =  1 ; 

pnOn  =  1 ; 

qbnfOn  =  1; 

wbnmOn  =  1 ; 

werrOn  =  1; 

berrOn  =  1; 

qerrOn  =  1; 
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%%  Set  Constants 
CONST. mu 

CONST . mu_moon  = 
CONST. mu_sun  = 
CONST. Re_ 

CONST. Rs 
CONST. J2 
CONST. J3 
CONST. J4 

CONST . SolarPress= 
CONST. SOLARSEC  = 
CONST. w_earth 
CONST. Cd 
CONST. Cr 
Reflect 

CONST . OmegaDot  = 
advance  for  sun-synch 


398 . 6004418el2;%mA3/sA2 
4 . 902 802 9535 97el2 ; %mA3/s A2 
1.327122E20; %mA3/sA2 
6.378137E6; %m 
1 . 4959787ell ; %m 
1. 0826266835 5E -3 ;% 
-2.53265648533E-6;% 

-1 . 61962159137E-6;% 

4 . 51e-6; %N/mA2 
806.81112382429; %TU 
-  [ 0 ; 0 ;  .000072  9211585530] ;%r/s 
2.5;% 

.6;% 


1. 991e-7;%rad/s 


earth  radius 
solar  radius 
J2  term 
J3  term 
J4  term 

solar  wind  pressure 

earth  rotation 
Coefficient  of  Drag 
Coefficient  of 

ascending  node 


%%  Set  Orbital  Elements 

%Kep  elements  meters  and  radians  (a, e, i, W, w, n) 


h_p  =  500e3;%m 

ha  =  500e3;%m 


altitude  at  perigee 
altitude  at  apogee 


RAAN  =  0 ; %rad 

w  =  0;%rad 

TAo  =  0 ; %rad 

Rp  =  CONST . Re +h_p; %m 

Ra  =  CONST . Re+h_a; %m 

e  =  (Ra-Rp) / (Ra+Rp) ; % (m/m) 

a  =  (Ra+Rp) /2;%m 

ho  =  sqrt (a*CONST .mu* ( l-eA2 ) ) ; %my2/s 

momentum 


Right  Ascention 
argument  of  perigee 
true  anomaly 
radius  of  perigee 
radius  of  apogee 
eccentricity 
semi -major  axis 
initial  angular 


P  =  2*pi* (aA3/CONST .mu) A . 5; %sec  Orbit  Period 

i_sunsynch  =  acosd ( (CONST . OmegaDot* ( l-eA2 ) A2 *aA ( 7 /2 )).. . 

/ ( -3/ 2 * sqrt (CONST . mu) * CONST . J2 * CONST . ReA2 ) ) ; %eqn  4.47  from 

Curtis 

i  =  i_sunsynch*D2R; %deg  (rad)  orbit  inclination 


[Ro,Vo]  =  sv_from_coe (CONST . mu, [ho  e  RAAN  i  w  TAo]);%  initial 
orbital  state  vector 


%%  Set  ICs 

w_BNo  =  [ 0 ; 2 *pi/P; 0 ] ; %rad  initial  body  rates 

w~ON  =  [  0 ; 2 *pi/P; 0 ] ; %rad 

rand ( ' seed ’ , 2 )  ; 
randn ( ' seed 1  ,  2 )  ; 
seedarw=l ; 
seedrrw=2 ; 
seedst=3 ; 
seedmag=4 ; 


108 


%  Sensor  parameters 
%  Gyro 

GYRO_Bias  =  ( 3*randn ( 3 , 1 ) ) *pi/180 ;  %  +  3  deg/sec 

N_ARW  =  (0 . 029) *pi/180; 

K_RRW  =  (0 . 0002) *pi/180; 

ARW  =  N_ARWA2;  %  angular  white  noise  Variance 

RRW  =  K_RRWA2 / 3 ;  %  bias  variance 

Gg  =  eye (3)  .* (-0.01  +  0. 02*rand(3) )  +... 

(ones (3,3) -eye (3)) .*(-0.0006+0. 0012*rand (3) ) ;  %percent 

%  Magnetometer 
sigMag  =  1 . 2  5  e - 7 ; 

Gm  =  eye (3) .* (-0.02+0. 04*rand(3) )  +... 

(ones (3,3) -eye (3)  )  .* (-0.0028  +  0. 0056*rand (3) ) ;  %percent 

%  Sun  Sensor 

51  =  [0  45  0] 1 *pi/180 ; 

52  =  [45  0  0] 1 *pi/180 ; 

SS_nl  =  [1  0  0]; 

SS_n2  =  [1  0  0]; 

FOV  =  0.7; 
sigSS  =  0.1; 

J  =  Bessel (sigSS/2, FOV) . *pi/180; 

%  Star  Tracker 

sigST  =  70  /3  /60  /60*pi/180;  %arcsec  to  rad  (3sig) 

%  Kalman  Filter 
dt  =  0.05; 
t_ekf  =  dt; 
sig(l)  =  sqrt (ARW) ; 
sig(2)  =  sqrt (RRW) ; 
sig(3)  =  sigST; 
sig(4)  =  sigSS*pi/180 ; 
sig(5)  =  sigMag; 

Ref erenceOmega  =  w_ON; 

[qBOo]  =  Euler_to_Quaternion ( InitialEuler ) ; 

[Ref erenceQuaternion]  =  Euler_to_Quaternion (Ref erenceEuler ) ; 

qBNo  =  qBOo; 

%%  Run  Simulation 
[Spacecraf t ] =  SCproperties ; 

J_Matrix  =  Spacecraf t . MOI ; 

[density_table]  =  GetDensity; 

RunTime  =  4000; %sec 


%sec  (20  Hz)  model  speed 
%sec  (100  Hz)  ekf  speed 
%rad/HzA (1/2) ,  ARW 
%rad/secA (3/2) ,  RRW 
%rad.  Star  Tracker  Error 
%rad.  Sun  Sensor  Error 
%tesla,  magnetometer  error 
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tic 

sim ( 1 ADS_Space craft Sim 1 , RunTime)  ; 

Total_Model_time  =  toe 

factor  =  RunTime/Total_Model_time 

DisturbanceTorques . Tgg  =  Tgg; 
DisturbanceTorques . Taero  =  Taero; 
DisturbanceTorques . Tsolar  =  Tsolar 

SensorMeasurements . ST  =  q_BNm; 
SensorMeasurements . Gyro  =  w_BNm; 
SensorMeasurements . bias  =  bias; 
SensorMeasurements . SSI  =  ssl; 
SensorMeasurements . SS2  =  ss2; 
SensorMeasurements . Mag  =  Bm; 

FilterEst.Q  =  q_BNf; 

FilterEst . Gyro  =  w_BNf; 

FilterEst . bias  =  bias  f; 


PlotUKFError.m 


bias_e_ul=squeeze  (bias_e_ul)  ; 

[m,n]=size (bias_e_ul ) ; 
if  m<n 

bias_e_ul=bias_e_ul ' ; 

end 

for  i=l : 3 

figure  ( 1 ) 
subplot (3, 1 , i) 

plot ( SimTimel , p_BNe_ul ( : , i ) ) 
hold  on 

plot (SimTimel , 3*sqrt ( Pdiag_ul ( : , i ) ) ,  '  -  .  r  '  ) 

plot (SimTimel , -3*sqrt ( Pdiag_ul ( : , i ) ) ,  '  -  .  r  '  ) 

grid  on 

xlim ( [0  4000] ) 

ylim ( [ - . 2  .2]) 

xlabel ( 'Time  (Sec) ' ) 

label= [ ' \deltap '  num2str(i)]; 

ylabel (label) 


end 

subplot (3,1,2) 
ylim ( [ - . 05  .05]) 
subplot  (3,1,1) 

title (' Generalized  Rodriquez  Parameter  Error  for  UKF ' ) 

for  i=4 : 6 

figure  ( 2 ) 
subplot  (3,1, i  —  3 ) 

plot (SimTimel , bias_e_ul ( : , i  —  3 ) ) ; 
hold  on 

plot (SimTimel , 3*sqrt (Pdiag_ul ( : , i) ) ,  ' - . r '  )  ; 

plot (SimTimel , -3*sqrt ( Pdiag_ul ( : , i ) ) , ' - . r ' ) ; 

grid  on 

xlim ( [0  4000] ) 

ylim ( [ -5E-4  5E-4] ) 

xlabel ( ' Time  (Sec) ' ) 

label= [ ' \delta  \beta  '  num2str ( i-3 ) ] ; 
ylabel (label) 

end 

subplot (3,1,1) 

title  ('Bias  Error  for  UKF') 


bias_e_el=squeeze (bias_e_el) ; 
[m, n] =size (bias_e_el ) ; 
if  m<n 


in 


bias_e_el=bias_e_el ' ; 

end 

for  i=l:3 

figure  ( 3 ) 
subplot (3, 1, i) 

plot ( SimTimel , q_BNe_el ( : , i ) ) 
hold  on 

plot (SimTimel , 3/2*sqrt (Pdiagl ( :  ,  i)  )  ,  '  -  .  r  '  )  ; 

plot (SimTimel , -3/2*sqrt (Pdiagl (  :  ,  i)  )  ,  1  - . r 1 ) 

grid  on 

xlim ( [0  4000] ) 

ylim ( [ - . 15  .15]) 

ylim ( [ - . 2  .2]) 

xlabel ( 1  Time  (Sec)  ' ) 

label= [ ' \deltaq '  num2str  (i) ] ; 

ylabel (label) 


end 

subplot (3,1,1) 

title (' Quaternion  Error  for  EKF') 
subplot (3,1,2) 
ylim ( [ - . 02  .02]) 

for  i=4 : 6 

figure  ( 4 ) 
subplot  ( 3 , 1 , i  —  3 ) 
plot (SimTimel, bias_e_el (:,i  —  3) ) 
hold  on 

plot (SimTimel , 3*sqrt (Pdiagl ( : , i) ) , '  -  .  r  '  )  ; 

plot (SimTimel , -3*sqrt (Pdiagl ( : , i) )  ,  '  -  .  r  '  )  ; 

grid  on 

xlim ( [0  4000] ) 

ylim ( [-5E-4  5E-4] ) 

xlabel ( 'Time  (Sec) ' ) 

label= [ ' \delta  \beta  '  num2str ( i-3 ) ] ; 
ylabel (label) 


end 

subplot (3,1,1) 

title  ('Bias  Error  for  EKF') 

clear  norm_p_u  norm_p_e  norm_bias_u  norm_bias_e 

norm_p_u=zeros (m,  1)  ; 

norm_p_e= zeros (m, 1) ; 

norm_bias_u= zeros (m,  1)  ; 

norm_bias_e=zeros (m,  1)  ; 

for  i=l : m; 

norm_p_u ( i , 1 ) =norm (p_BNe_ul ( i , : ) ) ; 
norm_p_e ( i , 1 ) =norm ( 2  *p_BNe_el ( i ,  : ) ) ; 
norm_bias_u ( i , 1 ) =norm (bias_e_ul ( i , : ) ) ; 
norm_bias_e ( i , 1 ) =norm (bias_e_el ( i , : ) ) ; 
end 

figure  ( 5 ) 
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semilogy ( SimTimel , norm_p_u) 
hold  on 

semilogy (SimTimel , norm_p_e, ' - . r ' ) 
grid  on 

title ( 1  Normalized  EKF  and  UKF  Attitude  Errors’) 
xlabel ( ' Time  (Sec) ’ ) 
ylabel (’ Attitude  Errors') 

legend (' Normalized  UKF  Generalized  Rodriguez  Parameter  Errors', 
'Normalized  EKF  Generalized  Rodriguez  Parameter  Errors') 
xlim ( [0  4000] ) 

figure ( 6) 

semilogy (SimTimel , norm_bias_u) 
hold  on 

semilogy (SimTimel , norm_bias_e,  ' - . r ' ) 
grid  on 

title (' Normalized  EKF  and  UKF  Bias  Errors') 

xlabel ( ' Time  (Sec) ' ) 

ylabel (' Normalized  \beta  Errors') 

legend (' Normalized  UKF  Bias  Errors',  'Normalized  EKF  Bias  Errors') 
ylim ( [ 0  .01]) 

xlim ( [0  4000] ) 


113 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


114 


LIST  OF  REFERENCES 


[1]  S.  Flagg,  R.  White,  and  R.  Ewart,  “Operationally  Responsive  Space  Specifications 
and  Standards:  An  Approach  to  Converging  with  the  Community,”  in  AIAA  Space 
2007  Conference  &  Exposition,  Long  Beach,  2007,  pp.  1-19. 

[2]  E.  J.  Lefferts,  F.  L.  Markley,  and  M.  D.  Shuster,  “Kalman  Filtering  for  Spacecraft 
Attitude  Estimation,”  Journal  of  Guidance,  Control  and  Dynamics,  vol.  15,  no.  5, 
1982,  pp.  417-429.  AIAA-82-0070. 

[3]  S.  J.  Julier  and  J.  K.  Uhlmann,  “Unscented  Filtering  and  Nonlinear  Estimation,” 
Proceedings  of  the  IEEE ,  vol.  92,  no.  3,  2004,  pp.  401-422. 

[4]  J.  L.  Crassidis  and  F.  L.  Markley,  “Unscented  Filtering  for  Spacecraft  Attitude 
Estimation,”  Journal  of  Guidance,  Control  and  Dynamics,  vol.  26,  no.  4,  2003,  pp. 
536-542. 

[5]  J.  Tuthill,  “Design  and  Simulation  of  a  Nano-Satellite  Attitude  Determination 
System,”  Naval  Postgraduate  School,  Monterey,  CA,  Master's  Thesis  December 
2009. 

[6]  L.  Stras,  D.  D.  Kekez,  et  al.  “The  Design  and  Operation  of  The  Canadian 
Advanced  Nanospace  experiment  (CanX-1),”  July  2004.  [Online],  Available: 
http://www.utias-sfl.net/nanosatellites/CanXl.  [Accessed  September  2010]. 

[7]  “AISsat-1  Monitoring  the  Shipping  Traffic  from  Space.”  NordicSpace  Online 
Journal,  April  2004.  [Online].  Available: 

http://www.nordicspace.net/PDF/NSA238.pdf.  [Accessed  September  2010]. 

[8]  M.  D.  Shuster  and  S.  D.  Oh,  “Three-Axis  Attitude  Determination  from  Vector 
Observations,”  Journal  of  Guidance,  Control  and  Dynamics,  vol.  4,  no.  1,  1981, 
pp.  70-77.  AIAA-8 1-4003. 

[9]  J.  L.  Crassidis  and  J.  L.  Junkins,  “ Optimal  Estimation  of  Dynamic  Systems ,” 

Goong  Chen  and  Thomas  J  Bridges,  Eds.  New  York,  USA:  Chapman  &  Hall/CRC, 
2004. 

[10]  M.  S.  Grewel  and  A.  A.  Andrews,  “ Kalman  Filtering:  Theory  and  Practice  Using 
MATLAB®”.  3rded.  John  Wiley  and  Sons,  2008. 

[11]  F.  Orderud,  “Comparison  of  Kalman  Filter  Estimation  Approaches  for  State  Space 
Models  with  Nonlinear  Measurements.,”  Sem  Scelands  vei,  vol.  7491,  pp.  7-9, 
March  2006. 


115 


[12]  S.  J.  Julier,  J.  K.  Uhlmann,  and  H.  F.  Durrant- Whyte.  “A  New  Approach  for 
Filtering  Nonlinear  Systems.”  In  Proceedings  of  the  1995  American  Control 
Conference,  1628-1632,  1995. 

[13]  M.  C.  VanDyke,  J.  L.  Schwartz,  and  C.  D.  Hall.  “Unscented  Kalman  Filtering  for 
Spacecraft  Attitude  State  and  Parameter  Estimation.”  In  Proceedings  of  the  2004 
Space  Flight  Mechanics  Meeting ,  AAS-04-1 15,  2004. 

[14]  A.  J.  Blocker,  “Tinyscope:  The  Feasibility  of  a  Tactically  Useful,  Three-Axis 
Stabilized,  Earth-Imaging  Nano-Satellite,”  Naval  Postgraduate  School,  Monterey, 
CA,  Master's  Thesis  December  2008. 

[15]  B.  Wei,  “Space  Vehicle  Dynamics  and  Control.”  AIAA  Education  Series,  1998. 

[16]  L.  F.  Markley,  “Multiplicative  vs.  Additive  Filtering  for  Spacecraft  Attitude 
Determination.”  July  2004.  [Online],  Available: 

http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov.  [Accessed  June  3,  2010]. 

[17]  S.  J.  Julier  and  J.  K.  Uhlmann,  “A  New  Extension  of  the  Kalman  Filter  to 
Nonlinear  Systems,”  In  Proceedings  of  the  SPIE  AeroSense  International 
Symposium  an  Aerospace/Defense  Sensing.  Simulations  and  Controls,  (Orlando, 
Florida),  April  20-25,  1997. 

[18]  H.  Curtis,  “Orbital  Mechanics  for  Engineering  Students,  2nd  ed.”  Elsevier’s 
Butterworth-Heinemann  2005. 

[19]  D.  Vallado,  “ Fundamentals  of  Astrodynamics  and  Applications,  2nd  ed.  ”  El 
Segundo:  Microcosm  Press,  2004. 

[20]  R.C.  Olsen,  “Introduction  to  the  Space  Environment.  ”  Monterey:  Naval 
Postgraduate  School,  2005. 

[21]  W.  Larson  and  J.  R.  Wertz,  “ Space  Mission  Analysis  and  Design,  3rd  ed.”  El 
Segundo:  Microcosm  Press,  2005. 

[22]  E.  Kraft,  “A  Quaternion-based  Unscented  Kalman  Filter  for  Orientation  Tracking,” 
Proceedings  of  the  6th  International  Conference  on  Information  Fusion,  Cairns, 
Australia,  2003,  pp.  47-54. 

[23]  C.  C.  Litton,  “Tinyscope:  The  Feasibility  of  a  Tactically  Useful  Earth-Imaging 
Nanosatellite  and  a  Preliminary  Design  of  the  Optical  Payload,”  Naval 
Postgraduate  School,  Monterey,  CA,  Master's  Thesis  January  2009. 


116 


INITIAL  DISTRIBUTION  LIST 


1 .  Defense  Technical  Information  Center 
Ft.  Belvoir,  Virginia 

2.  Dudley  Knox  Library 
Naval  Postgraduate  School 
Monterey,  California 


117 


